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I.  Introduction 


Predicting  the  behavior  of  combustion  waves  in  mixtures  of  gas  and  reactive  solid 
particles  is  an  important  and  partially  unsolved  problem.  Practical  applications  include  the 
burning  of  damaged,  granulated  solid  rocket  propellants,  detonation  of  granular  explosives, 
burning  of  coal  dust,  and  explosion  of  dust-air  mixtures.  Understanding  these  combustion 
processes  could  lead  to  more  accurate  design  criteria  for  rockets,  new  tailored  explosives, 
and  improved  safety  criteria  for  environments  where  dust  explosions  are  a  hazard. 

One  way  to  gain  understanding  is  to  model  these  processes.  A  class  of  models  which 
has  the  potential  to  describe  these  processes  has  been  developed  from  two-phase  continuum 
mixture  theory.  These  models  describe  each  phase  as  a  continuum;  distinct  equations  for 
the  mass,  momentum,  and  energy,  and  constitutive  equations  for  both  phases  are  written. 
The  two  phases  are  coupled  through  terms  representing  the  transfer  of  mass,  momentum, 
and  energy  from  one  phase  to  another.  Models  of  these  phase  interaction  processes  are 
determined  from  experiments.  In  the  models  the  phhse  interaction  terms  are  constructed 
such  that  global  conservation  of  mass,  momentum,  and  energy  is  maintained.  Regardless 
of  the  particular  form  of  the  two-phase  equations,  the  idea  of  global  conservation  is  a 
criterion  which  must  be  enforced. 

Unsteady  two-phase  models  have  been  widely  used  to  study  the  problem  of 
deflagration-to-detonation  transition  (DDT)  in  granulated  solid  propellants  [1-21],  which 
has  been  observed  experimentally  [22-24].  Similar  unsteady  models  are  used  to  study 
transient  combustion  in  porous  media  [25-30].  By  concentrating  on  unsteady  solutions, 
many  simple  results  available  from  the  less-complicated  two-phase  steady  theory  have  been 
overlooked.  These  results  are  found  by  solving  the  ordinary  differential  equations  which 
define  the  steady  two-phase  detonation  equilibrium  end  states  and  reaction  zone  structure. 
None  of  the  previous  steady  two-phase  studies  [31-39]  has  adequately  described  the 
admissible  end  states  and  structure  of  a  two-phase  detonation.  Only  when  steady 
detonation  solutions  are  understood  will  it  be  possible  to  fully  comprehend  the  implications 
of  unsteady  two-phase  detonation  theory. 

A  sketch  of  an  envisioned  two-phase  steady  detonation  structure  is  shown  in  Figure 
1.1.  In  this  study  the  term  "structure"  refers  to  the  spatial  details  of  the  detonation  wave. 
Such  details  include  the  reaction  zone  length  and  the  variation  of  pressure,  temperature,  etc. 
within  the  reaction  zone. 
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Figure  1.1  Hypothesized  Two-Phase  Steady  Detonation 

Drawing  on  the  results  of  one-phase  detonation  theory,  it  is  hypothesized  that  a  two-phase 
detonation  consists  of  a  chemicd  reaction  induced  by  a  shock  wave  propagating  into  a 
mixture  of  reactive  particles  and  inert  gas.  At  the  end  of  the  reaction  zone  the  particles  are 
completely  consumed;  only  inert  gas  remains.  Important  questions  concerning  such  a 
detonation  exist,  for  example, 

1)  What  is  the  speed  of  propagation  of  an  unsupported  two-phase 
detonation? 

2)  What  are  the  potential  two-phase  detonation  end  states? 

3)  What  is  the  structure  of  the  two-phase  reaction  zone? 

4)  What  is  the  nature  of  a  shock  wave  in  a  two-phase  material? 

5)  How  is  two-phase  detonation  theory  related  to  and  different  from  one- 
phase  detonation  theory? 

It  is  the  goal  of  this  work  to  use  steady  state  analysis  to  answer  these  and  other  questions. 

The  steady  equations  are  best  studied  using  standard  phase  space  techniques.  In  this 
work  such  techniques  are  used  to  study  a  general  two-phase  detonation  model.  In  so 
doing,  two-phase  steady  detonations  have  been  studied  in  the  same  context  as  the  extensive 
one-phase  steady  theory  [40]. 

An  outline  of  the  two-phase  detonation  analysis  of  this  work  is  now  given.  The 
unsteady  model  is  first  presented.  Then  the  steady  dimensionless  form  of  this  model  is 
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shown,  and  a  description  is  given  of  how  the  problem  of  determining  two-phase  detonation 
structure  can  be  reduced  to  solving  four  coupled  ordinary  differential  equations.  In  certain 
limits,  two  of  these  equations  may  be  integrated,  and  the  detonation  structure  problem  is 
reduced  to  solving  two  ordinary  differential  equations  in  two  unknowns.  In  these  limits  the 
detonation  structure  has  a  clear  geometrical  interpretation  in  the  two-dimensional  phase 
plane.  Both  two  and  four  equation  models  are  then  used  to  predict  examples  of  acceptable 
reaction  zone  structure  and  unacceptable,  non-physical  solutions.  Parametric  conditions  are 
obtained  for  the  existence  of  a  steady,  one-dimensional,  two-phase  detonation. 

Two  methods  are  used  to  restrict  the  available  solutions:  algebraic  end  state  analysis 
and  reaction  zone  structure  analysis.  Algebraic  analysis  of  the  equilibrium  end  states, 
described  in  detail  in  Ref.  41,  identifies  a  minimum  wave  speed  necessary  for  a  steady 
solution.  This  wave  speed  is  analogous  to  the  well-known  one-phase  Chapman-Jouguet 
(CJ)  wave  speed.  As  in  one-phase  theory,  the  two-phase  CJ  wave  speed  is  identified  as 
the  unique  wave  speed  of  an  unsupported  two-phase  detonation.  The  available  solutions 
are  further  restricted  by  considering  the  structure  of  the  two-phase  detonation  wave.  In 
particular,  results  from  the  stmctural  analysis  show  that  below  a  critical  initial  solid  volume 
fraction,  no  steady  two-phase  detonation  exists. 

The  behavior  of  integral  curves  near  singular  points  identified  by  this  analysis  is 
crucial  in  understanding  why  structural  analysis  limits  the  class  of  available  detonation 
solutions.  Analysis  of  two-phase  equations  near  singularities  has  not  been  emphasized  in 
two-phase  detonation  theory  or  two-phase  theory  in  general.  This  is  argued  by  Bilicki,  et 
al.  [42]  who  write  in  a  recent  article  concerning  steady  two-phase  flow, 

...the  theory  of  singular  points  of  systems  of  coupled,  ordinary  nonlinear 
differential  equations-still  largely  unexploited  in  this  field-is  essential  for  clarity, 
for  the  proper  management  ^computer  codes  and  for  the  understanding  of  the 
phenomenon  of  choking  as  predicted  by  the  adopted  mathematical  model,  an 
impossible  task  when  only  numerical  procedures  are  used. 

The  kingpin  of  the  analysis  is  the  identification  of  the  singular  points  of  the 
basic  system  of  equations  and  of  the  solution  patterns  that  they  imply.  Such  an 
analysis  serves  two  purposes.  First,  it  gives  the  analyst  the  ability  to  understand 
the  physical  characteristics  of  a  class  of  flows  without  the  need  to  produce 
complete  solutions.  Secondly,  it  gives  valuable  indications  as  to  how  to 
supplement  computer  codes  because  practically  all  numerical  methods  of  solution 
become  inadequate  in  the  neighborhood  of  the  singular  points  and  are 
constitutionally  incapable  of  locating  them  in  the  first  place,  which  leads  to 
numerical  difficulties  and  incorrect  interpretations.  This  has  to  do  with  the  fact 
that  the  set  of  algebraic  equations,  which  the  computer  code  must  solve  at  each 
step,  becomes  either  impossible  or  indeterminate... and  no  longer  solves  the 
coupled  differential  equations  of  the  model. 
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The  analysis  presented  here  identifies  two  types  of  singular  points,  explained  in  detail 
below,  which  exist  in  most  two-phase  particle  burning  models  based  on  continuum  mixture 
theory.  Near  a  singularity  there  is  a  zero  in  the  denominator  of  the  forcing  functions  of  the 
governing  differential  equations.  The  consequences  of  these  singularities  are  not 
straightforward  and  must  be  analyzed  in  detail. 

One  type  of  singularity  occurs  at  the  point  of  complete  reaction.  The  complete  reaction 
singularity  arises  in  most  particle- burning  two-phase  detonation  models  because  the 
interphase  transport  terms  used  in  the  mass,  momentum,  and  energy  equations  typically 
have  a  1/r  dependence  where  r  is  particle  radius.  When  the  particle  radius  approaches  zero, 
a  singularity  exists.  It  is  an  open  question  as  to  whether  this  singularity  gives  rise  to 
infinite  gradients,  infinite  property  values,  or  whether  there  is  a  balancing  zero  in  the 
numerator  to  counteract  the  singularity.  No  work  in  the  current  two-phase  detonation 
literature  adequately  addresses  this  issue.  The  results  presented  in  this  work  account  for 
the  complete  reaction  singularity. 

Another  type  of  singularity  occurs  when  the  velocity  of  either  phase  relative  to  the 
wave  front  is  locally  sonic.  In  this  work  the  term  "sonic"  is  taken  to  mean  that  the  velocity 
of  an  individual  phase  relative  to  the  steady  wave  is  equal  to  the  local  sound  speed  of  that 
particular  phase  as  predicted  by  the  state  equation  for  that  particular  phase.  The  term 
"sonic"  in  this  work  does  not  in  any  way  refer  to  a  mixture  sound  speed,  nor  is  the  idea  of 
a  mixture  sound  speed  incorporated  into  any  of  the  arguments  developed  in  this  work. 

The  sonic  singularity  arises  naturally  from  the  differential  equations  and  has  been 
extensively  studied  for  one-phase  systems.  Here  for  the  first  time  the  importance  of  sonic 
conditions  in  two-phase  detonation  theory  is  shown:  in  general  if  a  solid  sonic  condition  is 
reached  within  the  detonation  structure,  a  physically  acceptable  steady  two-phase 
detonation  cannot  exist.  If  a  solid  sonic  condition  is  reached,  it  is  predicted  that  all  physical 
variables  are  double-valued  functions  of  position;  for  instance  at  any  point  in  the  wave 
structure  two  distinct  gas  densities,  solid  temperatures,  etc.  are  predicted.  This  condition  is 
obviously  not  physical.  Furthermore,  when  such  a  condition  is  reached  the  solution  does 
not  reach  an  equilibrium  point;  thus,  no  steady  solution  is  predicted.  This  alone  is  a 
sufficient  reason  to  reject  solutions  contain  a  solid  sonic  condition.  In  addition  to  the  solid 
phase  sonic  singularity,  imaginary  gas  phase  properties  are  predicted  if  the  solution 
includes  a  gas  phase  sonic  point  at  a  point  of  incomplete  reaction. 

The  influence  of  the  two-phase  shock  state  on  steady  detonation  stmcture  is  shown  in 
this  work.  The  shock  wave,  assumed  to  be  inert,  leaves  the  material  in  a  state  of  higher 
pressure,  temperature,  and  density  than  the  ambient  slate.  This  serves  to  initiate  chemical 
reaction  which  in  turn  releases  energy  to  drive  the  shock  wave.  In  this  work  mechanisms 
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which  define  the  structure  of  a  shock  wave  such  as  diffusive  heat  conduction  and 
momentum  transport  are  ignored.  It  is  assumed  that  the  length  scales  on  which  these 
processes  are  important  are  small  in  comparison  with  the  reaction  zone  length  scales.  By 
ignoring  tb:  diffusive  processes,  the  model  equations  become  hyperbolic,  and 
discontinuous  shocks  are  admitted  by  the  governing  equations. 

i  ne  shock  state  can  be  determined  by  an  algebraic  analysis.  Any  state  admitted  by  the 
shock  discontinuity  equations  can  serve  as  an  initial  condition  for  the  ordinary  differential 
equations  which  define  the  reaction  zone  structure.  It  is  shown  that  four  classes  of  initial 
conditions  are  admitted  for  a  given  wave  speed:  1)  gas  and  solid  at  ambient  conditions,  2) 
a  shocked  gas  and  shocked  solid,  3)  an  unshocked  gas  and  shocked  solid,  and  4)  a 
shocked  gas  and  unshocked  solid.  Any  of  these  initial  states  has  the  potential  to  initiate  a 
steady  two-phase  detonation.  Examples  are  found  of  the  fu'st  and  fourth  classes  of  two- 
phase  detonation  in  this  thesis.  Previous  work  in  two-phase  detonation  has  not  adequately 
shown  whether  the  gas  and  solid  are  shocked  or  unshocked. 

In  addition  to  two-phase  detonation  structure,  this  study  contains  a  discussion  of  inert 
compaction  waves  in  granular  materials.  This  discussion,  including  a  review  of 
compaction  wave  theory  and  experiments,  is  contained  in  Chapter  4  and  is  not  germane  to 
the  subject  of  steady  two-phase  detonations.  The  results  are  predicted  by  the  same 
equations  used  to  predict  two-phase  detonations  in  the  limit  of  no  chemical  reaction  and  a 
negligible  gas  phase.  In  Chapter  4  analysis  is  presented  to  describe  the  wave  motion  which 
results  when  a  constant  velocity  piston  strikes  a  granular  material. 

A  sketch  of  an  envisioned  compaction  wave  is  shown  in  Figure  1 .2. 


Compaction  Wave  Speed  =  D 


^  I  Undisturbed  Region 

Compaction 

Zone 

Figure  1.2  Sketch  of  Compaction  Wave  in  Granular  Material 
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A  compaction  wave  is  thought  to  be  an  important  event  in  the  process  in  the  transition 
from  deflagration  to  detonation  (DDT).  It  is  thought  that  a  compaction  process  in  which  the 
granular  material  rearranges  can  give  rise  to  local  hot  spots  which  could  induce  a  detonation 
in  the  reactive  material. 

Much  as  for  two-phase  steady  detonation  analysis,  the  compaction  wave  analysis 
identifies  equilibrium  end  states  and  compaction  zone  structure.  It  is  shown  that  the 
problem  of  determining  compaction  zone  structure  can  be  reduced  to  solving  one  ordinary 
differential  equation  for  one  unknown,  solid  volume  fraction.  The  results  show  a 
continuous  dependence  of  compaction  wave  structure  with  supporting  piston  velocity; 
depending  on  the  piston  velocity,  two  broad  classes  of  compaction  zone  structure  exist.  At 
low  piston  velocities  the  compaction  wave  travels  at  speeds  less  than  the  ambient  solid 
sound  speed.  Such  waves  are  called  subsonic  compaction  waves.  The  structure  is 
characterized  by  a  smooth  rise  in  pressure  from  the  ambient  to  a  higher  pressure  equal  to 
the  static  pore  collapse  stress  level.  Subsonic  compaction  waves  have  been  observed  in 
experiment  [43,  44]  and  predicted  by  Baer  [45]  and  Powers,  Stewart,  and  Krier  [46]. 
Above  a  critical  piston  velocity  the  compaction  wave  travels  at  speeds  greater  than  the 
ambient  solid  sound  speed.  A  discontinuous  shock  wave  leads  a  relaxation  zone  where  the 
pressure  adjusts  to  its  equilibrium  static  pore  collapse  value.  Such  waves  are  called 
supersonic  compaction  waves.  Supersonic  compaction  waves  with  leading  shocks  have 
not  yet  been  observed  nor  predicted  in  previous  studies. 

The  plan  of  this  thesis  is  to  first  review  the  relevant  literature  in  Chapter  2.  The 
unsteady  model  is  presented  in  Chapter  3.  Steady  inert  compaction  waves  predicted  by  this 
model  are  discussed  in  Chapter  4  which  is  followed  by  a  discussion  of  two-phase 
detonation  equilibrium  end  states  and  structure  in  Chapter  5.  Conclusions  and 
recommendations  are  given  in  Chapter  6.  Appendix  A  discusses  the  method  of 
characteristics,  and  lists  the  characteristic  directions  and  equations  for  one-dimensional, 
unsteady,  two-phase  reactive  flow.  Appendix  B  has  a  detailed  discussion  of  state  relations 
and  demonstrates  that  the  thermal  and  caloric  state  equations  used  in  this  study  are 
compatible.  Appendix  C  compares  the  momentum  and  energy  equations  of  this  study  to 
other  common  forms  of  these  equations  and  defends  the  choices  made  for  this  study.  Two- 
phase  CJ  deflagration  conditions  are  considered  in  Appendix  D.  Appendix  E  contains  a 
detailed  description  of  how  to  reduce  the  model  to  the  simple  two-equation  model  presented 
in  Chapter  5.  Appendix  F  gives  a  derivation  of  the  number  conservation  equation.  This 
equation  holds  that  in  the  two-phase  flow  field,  the  number  density  of  particles  does  not 
change. 
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n.  Review  of  Two-Phase  Detonation  Theory 


This  chapter  will  briefly  describe  the  literature  which  is  relevant  to  the  field  of  two- 
phase  steady  detonation  theory.  This  includes  works  on  the  fundamentals  of  two-phase 
continuum  mixture  theory,  basic  one-phase  detonation  theory,  and  applications  of  these 
theories  to  combustion  in  porous  media.  As  this  thesis  is  primarily  concerned  with  the 
details  of  modeling  two-phase  detonations  using  existing  models  and  not  with  the 
experiments  which  provide  the  basis  for  these  models,  the  experimental  literature  regarding 
two-phase  detonations  will  not  be  reviewed.  The  interested  reader  is  referred  to  Butler’s 
thesis  [47]  for  a  thorough  description.  A  review  of  compaction  wave  theory  is  found  in 
Chapter  4. 

The  theory  of  two-phase  flow  is  still  under  development,  and  there  are  many  issues 
which  remain  unresolved.  Drew  [48]  considers  some  of  these  issues  in  a  recent  review 
article.  However,  one  need  only  look  at  the  wide  disparity  in  the  forms  of  two-phase 
model  equations  expressed  by  various  researchers  to  realize  that  the  particular  form  of  the 
equations  is  a  matter  of  dispute.  Thus  in  constructing  a  model,  one  looks  for  the  most 
basic  principles  to  use  as  a  guide.  In  his  description  of  two-phase  theory  from  a  continuum 
mechanics  perspective,  Truesdell  [49]  describes  three  metaphysical  principles  which  can  be 
used  as  a  guide.  They  are: 

1.  All  properties  of  the  mixture  must  be  mathematical  consequences  of 
properties  of  the  constituents. 

2.  So  as  to  describe  the  motion  of  a  constituent,  we  may  in  imagination 
isolate  it  from  the  rest  of  the  mixture,  provided  we  allow  properly  for 
the  actions  of  the  other  constituents  upon  it. 

3.  The  motion  of  the  mixture  is  governed  by  the  same  equations  as  is  a 
single  body. 

Two-phase  theory  as  applied  to  combustion  in  granular  materials  has  been  developed 
primarily  through  the  work  of  Krier  and  co-workers  [1,  7,  16,  17,  18,  20,  21],  Kuo, 
Summerfield,  ana  co-workers  [29,  30,  37,  38],  and  more  recently  by  Nunziato,  Baer,  and 
co-workers  [2,  3,  4,  11,  13].  In  addition,  Nigmatulin's  book  [50],  available  in  Russian 
■id  not  reviewed  by  this  author,  is  widely  referred  to  in  the  Russian  literature  as  a  source 
for  the  governing  equations  of  two-phase  reactive  flow.  As  opposed  to  the  work  of  Kuo, 
et  al.,  who  consider  only  deflagrations,  the  work  of  Krier,  et  al.,  and  Nunziato,  et  al.,  has 
been  applied  to  the  detonation  of  granular  propellants  and  explosives.  The  extreme 


conditions  of  a  detonation  (gas  pressures  are  of  the  order  of  10  GPa)  force  the  adoption  of 
a  fully  compressible  solid  phase  state  equation  and  a  non-ideal  gas  phase  state  equation. 
These  models  use  constitutive  theory  specially  developed  to  describe  the  pore  collapse 
which  can  be  associated  with  detonations  in  these  materials.  There  are  only  minor 
differences  in  the  Krier  and  Nunziato  model  formulations;  these  are  considered  in  detail  in 
Chapter  3  and  Appendix  C. 

To  understand  two-phase  detonation  theory,  it  is  necessary  to  be  familiar  with  some  of 
the  results  of  one-phase  detonation  theory.  The  best  summary  of  these  results  is  given  in 
Pickett  and  Davis's  book  [40].  The  one-phase  concept  most  relevant  to  two-phase  theory 
is  that  of  a  steady  Zeldovich,  von  Neumann,  Doering  (ZND)  detonation  which  terminates  at 
a  CJ  point.  The  ZND  theory  is  named  for  its  developers  who  independently  described  the 
theory  in  the  1940's.  The  CJ  analysis  describes  the  equilibrium  end  states  for  a  one-phase 
detonation,  and  the  ZND  analysis  describes  the  structure  of  the  detonation  reaction  zone. 

Much  of  one-phase  detonation  theroy  can  be  understood  by  considering  the 
equilibrium  end  states.  The  (3J  point  is  an  equilibrium  end  state  at  which  the  gas  velocity  is 
sonic  with  respect  to  the  wave  front.  Since  this  point  is  sonic,  the  theory  predicts  that  any 
trailing  rarefaction  wave  is  unable  to  catch  and  disturb  the  steady  wave.  There  is  only  one 
detonation  wave  speed  which  leaves  the  material  in  a  CJ  state.  The  equilibrium  end  state 
analysis  of  one-dimensional  theory  hypothesizes  that  this  wave  speed  is  the  unique  speed 
of  propagation  for  an  unsupported  detonation  wave.  There  are  no  equilibrium  end  states 
for  wave  speeds  less  than  the  CJ  wave  speed.  For  wave  speeds  greater  than  the  CJ  wave 
speed,  two  equilibrium  states  are  predicted.  They  are  classified  on  the  basis  of  the 
equilbrium  end  state  pressure:  the  solution  which  terminates  at  the  higher  pressure  is  called 
a  strong  solution,  the  other  solution  is  called  the  weak  solution.  The  strong  end  state  is  a 
subsonic  state,  and  thus  the  strong  detonation  is  susceptable  to  degradation  from  trailing 
rarefactions.  To  achieve  a  strong  detonation,  the  theory  predicts  a  supporting  piston  is 
necessary  so  that  no  rarefactions  will  exist.  The  weak  end  state  is  a  supersonic  state  and 
thus  does  not  require  any  piston  support  and  is  not  ruled  out  by  simple  equilibrium  end 
state  analysis. 

ZND  theory  considers  the  structure  of  a  detonation  wave  which  links  the  initial  state  to 
the  equilibrium  end  state.  A  ZND  detonation  is  described  by  an  inert  shock  wave 
propagating  into  a  reactive  material.  The  shock  wave  leaves  the  material  in  a  locally 
subsonic,  high  temperature  state.  The  high  temperature  initiates  an  exothermic  chemical 
reaction.  Energy  released  by  this  chemical  reaction  is  predicted  to  drive  the  detonation 
wave.  For  wave  speeds  greater  than  the  CJ  speed,  the  solution  terminates  at  the  strong 
point,  a  subsonic  state  which  requires  piston  support  to  remain  steady.  For  a  CJ  wave 


speed  the  solution  terminates  at  a  sonic  point  and  thus  is  able  to  propagate  without  piston 
support.  The  simple  ZND  theroy  predicts  that  there  is  no  path  from  the  initial  shock  state  to 
the  weak  solution  point  and  thus  rules  out  a  weak  detonation  with  a  leading  shock  in  the 
structure.  Thus  simple  ZND  theory  predicts  that  the  wave  speed  for  an  unsupported 
detonation  is  the  O  wave  speed.  There  is,  however,  evidence,  described  in  detail  by 
Pickett  and  Davis,  that  weak  solutions  can  be  achieved.  In  general  the  weak  detonations 
described  by  Pickett  and  Davis  require  special  conditions  to  exist 

Pickett  and  Davis  describe  how  ZND  theory  can  be  placed  in  the  context  of  the  general 
theory  of  systems  of  ordinary  differential  equations.  Details  of  this  theory  can  be  found  in 
standard  texts  [51,  52].  The  theory  describes  how,  given  a  set  of  ordinary  differential 
equations,  solutions  link  an  initial  state  to  an  equilbrium  end  state  and  how  other  solutions 
do  not  have  equilibrium  states.  Equilibrium  states  are  defined  at  points  where  the  forcing 
functions  for  each  differential  equation  are  simultaneously  zero.  Whether  or  not  an 
equilibrium  state  is  reached  depends  on  the  particular  form  of  the  differential  equations.  A 
solution  which  does  not  reach  an  equilibrium  point  is  rejected  as  a  steady  solution  by 
definition. 

A  shortcoming  of  most  two-phase  detonation  studies  is  that  little  emphasis  has  been 
put  on  placing  two-phase  detonation  theojy  in  the  context  of  one-phase  detonation  theory 
and  the  more  general  ordinary  differential  equation  theory.  Most  work  in  two-phase 
detonation  theory  has  concentrated  solely  on  using  numerical  solution  of  the  unsteady 
equations  to  predict  the  two-phase  equivalent  of  a  CJ  detonation  [1,2].  A  primary  goal  of 
these  works  has  been  to  predict  the  deflagration-to-detonation  transition  (DDT)  zone  length 
rather  than  the  character  of  the  detonation  itself.  As  such,  there  has  been  little  discussion  of 
the  basic  characteristics  of  a  steady  two-phase  detonation.  In  these  studies  the  definition  of 
the  two-phase  CJ  state  is  unclear.  Reference  is  often  made  to  the  one-phase  CJ  results  with 
the  assumption  that  the  one-phase  CJ  condition  naturally  must  also  apply  to  the  two-phase 
detonation  model.  Also  in  these  works  no  attempt  has  been  made  to  describe  conditions 
under  which  a  two-phase  strong  or  weak  detonation  can  exist.  Since  these  states  can  be 
predicted  by  one-phase  theory,  it  is  reasonable  to  suggest  that  two-phase  equivalents  may 
also  exist.  Detailed  descriptions  of  the  steady  two-phase  reaction  zone  structure  have  been 
generally  ignored. 

Studies  of  steady  two-phase  systems  will  now  be  considered.  Several  works  exist 
which  consider  the  relatively  low-speed,  low-pressure  deflagration  of  solid  particles. 
Among  these  are  the  works  of  Kuo,  et  al.  [37,  38],  Ermolaev,  et  al.  [35,  36],  and  Drew 
[31].  These  works  consider  the  particle  phase  to  be  incompressible  and  naturally  have  no 
discussion  of  shock  waves.  The  work  of  Krier  and  Mozafarrian  [34]  considered  a  reactive 
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wave  with  a  leading  shock  wave  in  the  gas  phase.  Detonation  structure  was  determined  by 
numerically  solving  the  steady  two-phase  ordinary  differential  equations.  This  work  is  of 
limited  value  because  of  the  assumption  of  an  incompressible  solid  phase.  This  assumption 
is  unrealistic  in  the  detonation  regime.  In  addition  they  did  not  establish  whether  the  model 
equations  of  their  work  are  hyperbolic,  leading  one  to  question  whether  their  model 
equations  are  well-posed  for  their  initial  value  problem. 

The  most  important  studies  of  steady  two-phase  detonation  are  those  of  Sharon  and 
Bankoff  [33]  and  Condiff  [32].  These  works  apply  two-phase  detonation  theory  to  vapor 
explosions  which  can  arise  from  the  rapid  mixing  of  a  hot  liquid  and  cold  vaporizable 
liquid.  Large  differences  in  the  features  of  the  problem  of  vapor  explosion  and  that  of 
detonation  of  solid  granular  explosive  prevent  an  extension  of  results  of  vapor  explosion  to 
detonations  in  granular  explosives  from  being  made.  Among  the  differences  are  that  in  a 
vapor  explosion  both  components  come  to  an  equilibrium  where  both  components  exist  in 
finite  quantities,  while  in  a  detonation  of  a  solid  propellant  the  solid  is  entirely  consumed. 
There  are  also  large  differences  in  the  functional  form  of  the  constitutive  equations. 
Nevertheless,  both  these  works  discuss  in  detail  many  features  of  two-phase  detonation 
theory  which  are  held  in  common  between  vapor  explosions  and  detonations  in  granular 
explosives.  More  importandy,  these  works  oudine  a  rational  approach  to  the  problem  of 
two-phase  detonation. 

Both  Sharon  and  Bankoff  and  Condiff  describe  a  two-phase  detonation  in  the  context 
of  one-phase  steady  detonation  theory.  That  is  they  describe  the  detonation  structure  as  a 
shock  jump  followed  by  a  relaxation  zone  whose  structure  is  determined  by  solving  a  set  of 
ordinary  differential  equations.  Both  works  describe  the  effective  two-phase  CJ  state. 
Sharon  and  Bankoff  argue  that  the  CJ  vapor  explosion  is  the  only  steady  solution  which 
can  exist.  They  also  provide  details  of  the  detonation  structure.  Condiff  argues  that  there 
is  a  larger  range  of  solutions  to  choose  from  and  that  a  phase  plane  analysis  is  necessary  to 
choose  which  solutions  can  be  accepted. 

An  issue  which  has  long  been  troublesome  for  two-phase  theory  is  that  of  whether  the 
equations  are  well-posed.  In  two-phase  detonation  theory,  only  two  models  have  been 
prop)osed  which  have  been  shown  to  be  well-posed  for  initial  value  problems:  the  model  of 
Baer  and  Nunziato  [2]  and  Powers,  Stewart,  and  Krier  [41] .  The  feature  of  these  models 
which  guarantees  that  they  are  well-posed  is  an  explicit  time-dependent  equation  which 
models  the  change  in  volume  fraction  in  a  granular  material.  When  models  without  such  an 
equation  are  examined,  it  is  found  that  there  are  regimes  in  which  imaginary  characteristics 
arc  present  [8,  28,  48,  53].  Such  models  are  not  in  general  well-posed  for  initial  value 


11 


problems;  because  of  this  any  solution  to  an  initial  value  problem  for  such  a  model  can  be 
shown  to  be  unstable  to  disturbances  of  any  ftequency. 

It  should  be  said  that  for  gas  phase  systems  that  the  ZND  assumption  of  one- 
dimensionality  has  been  shown  by  experiments  to  be  invalid  in  general.  However  the  ZND 
predictions  are  able  to  roughly  predict  spatially  averaged  gas  phase  properties  such  as  final 
pressure  and  wave  speed.  For  solids,  experimental  results  provide  little  evidence  regarding 
the  existence  of  multidimensional  detonation  structure.  Regardless  of  whether  or  not 
detonations  in  solids  are  one  or  multidimensional,  it  is  reasonable  to  consider  the  results  of 
one-dimensional  theory  before  proceeding  to  consider  more  complicated  multidimensional 
theories. 

Finally  an  issue  relevant  to  models  of  particle  burning  must  be  considered,  that  of  how 
an  expression  for  the  evolution  of  particle  radius  should  be  formulated.  In  two-phase 
models  of  granular  materials  the  particle  radius  is  a  required  variable  for  all  interphase 
transfer  terms  (reaction,  drag,  and  heat  transfer  are  known  empirically  as  functions  of 
particle  radius).  In  much  of  the  two-phase  granular  explosive  literature  there  is  confusion 
as  to  how  to  determine  the  particle  radius.  The  recent  work  of  Baer  and  Nunziato  [2] 
disregards  the  issue  by  not  giving  an  expression  for  particle  radius  evolution.  It  would 
seem  that  this  model  is  incomplete.  The  work  of  Krier  and  co-workers  provides  a  relation 
for  the  particle  radius  evolution  whose  physical  interpretation  is  unclear  (see  Appendix  F). 
A  rational  way  for  determining  particle  radius  is  found  by  considering  an  evolution 
equation  for  the  number  density  of  particles.  Several  modelers  do  write  explicit  equations 
for  number  density  evolution  [10, 14, 26, 31, 39].  Generally  these  studies  assume  that  the 
number  density  of  particles  is  conserved.  It  can  be  shown  that  with  such  an  equation  it  is 
possible  to  determine  a  clearly  understood  equation  for  the  evolution  of  particle  radius  (see 
Appendix  F). 
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m.  The  Unsteady  Two-Phase  Model 


A  two-phase  model  is  presented  which  is  a  slight  modification  of  the  model  first 
presented  in  Ref.  41.  It  is  similar  to  models  used  by  Butler  and  Krier  [1]  and  Baer  and 
Nurmato  [2].  Changes  of  two  types  have  been  made.  First  a  simplified  constitutive  theory 
has  been  adopted  in  order  to  make  the  equations  tractable.  The  trends  predicted  by  the 
simpler  constitutive  equations  are  similar  to  the  trends  of  Refs.  1  and  2.  A  second  more 
substantial  change  is  that  an  explicit  expression  for  particle  radius  evolution  has  been 
adopted.  No  counterpart  to  this  equation  is  found  in  either  Ref.  1  or  2.  For  the  proposed 
model  it  is  assumed  that  each  phase  is  a  continuum;  consequently,  partial  differential 
equations  resembling  one-phase  equations  are  written  to  describe  the  evolution  of  mass, 
momentum,  and  energy  in  each  constituent.  In  addition,  each  phase  is  described  by  a 
thermal  state  relation  and  a  corresponding  caloric  state  relation.  Constituent  one  is  assumed 
to  be  a  gas,  constiment  two,  a  solid. 

In  order  to  close  the  system,  a  dynamic  compaction  equation  similar  to  that  of  Ref.  2  is 
adopted.  Choosing  a  dynamic  compaction  equation  insures  that  the  characteristics  are  real; 
thus,  the  initial  value  problem  is  well-posed.  The  unsteady  two-phase  model  is  posed  in 
characteristic  form  in  Appendix  A.  The  dynamic  compaction  equation  states  that  the  solid 
volume  fraction  changes  in  response  to  1)  a  difference  between  the  solid  pressure  and  the 
sum  of  the  gas  pressure  and  intragranular  stress  and  2)  combustion.  Most  models  use 
empirical  data  to  model  the  intragranular  stress.  Here  for  simplicity  it  is  assumed  that  the 
intragranular  stress  is  a  linear  function  of  the  solid  volume  fraction.  This  function  is 
constructed  such  that  no  volume  fraction  change  due  to  pressure  differences  is  predicted  in 
the  initial  state.  It  is  emphasized  that  the  choices  made  for  the  closure  problem  and  for 
other  constitutive  relations  place  a  premium  on  simplicity  so  that  explicit  analytic 
calculations  can  be  made  whenever  possible.  At  the  same  time  the  model  adopted  here  is 
representative  of  a  wider  class  of  two-phase  detonation  models  currently  in  use. 

The  unsteady  equations  are 
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Here  the  subscript  "0"  denotes  the  undisturbed  condition,  "1"  denotes  the  gas  phase; 
"2,"  solid  phase;  p ,  density;  <(>,  volume  fraction;  u,  velocity;  r,  solid  particle  radius;  P, 
pressure;  m,  bum  index;  a,  bum  constant;  P,  drag  parameter,  e,  internal  energy;  h,  heat 
transfer  coefficient;  R,  gas  constant;  b,  co- volume  correction;  c,  sound  speed;  c^,  constant 
volume  specific  heat;  s,  non-ideal  solid  parameter;  Pc  compaction  viscosity;  Tait 
parameter,  and  q,  heat  of  reaction. 

Numerical  values  for  the  parameters  introduced  above,  representative  of  the  solid  high 
explosive  HMX,  are  listed  in  Table  1.  When  available,  references  are  listed  for  each  of  the 
parameters.  The  unreferenced  parameters  have  been  estimated  for  this  study.  The  initial 
gas  density  and  temperature  have  arbitrarily  been  chosen  to  be  10  kg/m^  and  300  K, 
respectively.  Drag  and  heat  transfer  parameters  have  been  chosen  to  roughly  match 
empirical  formulae  given  in  Ref.  13.  The  gas  constant  R  and  virial  coefficient  b  have  been 
chosen  such  that  predictions  of  CJ  detonation  states  match  the  CJ  detonation  states 
predicted  by  the  thermochemistry  code  TIGER  [54]  as  reported  in  Ref.  1.  The  solid 
parameters  s  and  Y2  have  been  chosen  such  that  solid  shock  and  compaction  wave 
predictions  match  experimental  shock  [55]  and  compaction  wave  data  [43,  44].  As 
reported  by  Baer  and  Nunziato  [2],  there  are  no  good  estimates  for  the  compaction 
viscosity  Ref.  2  chooses  a  value  for  compaction  viscosity  of  10^  kg/m  s.  To 
demonstrate  the  existence  of  a  two-phase  detonation,  it  was  necessary  in  this  study  to 
choose  a  higher  value,  10^  kg/m  s,  for  the  compaction  viscosity. 

Undisturbed  conditions  are  specified  as 


p  =  p 

^10 


p  =  p 
*^2  ^20 


(j)  =  <j) 

1  10 


u^  =  0, 


u 


2 


=  0 


T 


Undisturbed  conditions  for  other  variables  can  be  determined  by  using  the  algebraic 
relations  (3.9-15). 

Equations  (3.1,2)  describe  the  evolution  of  each  phase's  mass;  Equations  (3.3,4), 
momentum  evolution;  and  Equations  (3.5,6),  energy  evolution.  Homogeneous  mixture 
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Table  I 

Dimensional  Input  parameters 


a 

[1] 

[m/(sPa)] 

2.90  X  10-9 

Pio 

[kg/ m3] 

1.00  X  101 

m 

[1] 

1.00  X  lOO 

P 

[kg/(sm2)] 

1.00  X  10^ 

P20 

[1.2] 

[kg/ m3] 

1.90  X  103 

h 

[J/(sKm8^)] 

1.00  X  102 

Cvl 

[2] 

[J/(kgK)] 

2.40  X  103 

Cv2 

[1,2] 

[J/(kgK)] 

1.50  X  103 

R 

[J/(kgK)] 

8.50  X  102 

s 

[(m/s)2] 

8.98  X  106 

q 

[1] 

[J/kg] 

5.84  X  106 

ro 

[1.2] 

[m] 

1.00  X  10-4 

b 

[m3 /kg] 

1.10  X  10-3 

Y2 

5.00  X  lOO 

[kg /(ms)] 

1.00  X  106 

To 

[K] 

3.00  X  102 
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equations  are  formed  by  adding  Equations  (3.1)  and  (32),  (3.3)  and  (3.4),  and  (3.5)  and 
(3.6).  Thus  for  the  mixture,  conservation  of  mass,  momentum,  and  energy  is  maintained. 

The  forcing  functions,  inhomogeneities  in  Equations  (3.1-6),  model  inter-phase 
momenti:m,  energy,  and  mass  transfer.  Functional  forms  of  inter-phase  transfer  terms 
have  been  chosen  to  have  a  simple  form.  Figure  3.1  shows  a  comparison  of  the  above 
drag  model  and  the  empirical  model  used  by  Baer  [13],  which  is  dependent  on  Reynolds 
number  for  particle  radii  from  0  to  300  jxm.  The  Reynolds  number  has  been  found  to  lie  in 
the  range  0-1000  within  the  two-phase  detonation  reaction  zones  of  this  study.  Figure  3.1 
shows  that  the  functional  forms  of  the  two  relations  are  similar,  though  the  magnitudes 
vary  widely.  A  similar  comparison  is  made  in  Figure  3.2  between  the  simplified  inter¬ 
phase  heat  transfer  modelled  here  and  the  empirical  heat  transfer  model  used  by  Baer. 
Again,  the  functional  form  of  the  two  models  is  similar  and  wide  variation  exists  in  the 
magnitudes.  For  mass  transfer  a  well-known  empirical  relation  for  the  regression  of 
particle  radius  is  used.  It  is  observed  that  the  rate  of  change  of  particle  radius  is 
proportional  to  the  surrounding  pressure  raised  to  some  power.  The  right  sides  of  the  mass 
equations  (3.1-2)  are  formulated  to  incorporate  this  feature. 

By  combining  the  solid  mass  evolution  equation  (3.2)  with  the  number  conservation 
equation  (3.7),  an  explicit  equation  is  obtained  fcMr  particle  radius  evolution: 


dr  dr 

at  “2  ax 


2  ax 


(3.16) 


This  equation  demonstrates  that  following  a  particle,  the  particle  radius  may  change  in 
response  to  combustion,  embodied  in  the  empirically- based  term  -aPi"*,  and  density 
changes,  as  described  by  the  density  derivative  terms.  Many  two-phase  particle-burning 
detonation  models  do  not  explicitly  include  an  equation  for  the  evolution  of  particle  radius. 
In  these  models,  which  also  do  not  explicitly  enforce  number  conservation,  it  is  unclear 
what  physical  principles  are  used  to  determine  the  particle  radius.  For  a  detailed  derivation 
of  the  number  conservation  relation  (3.7)  and  Equation  (3.16)  see  Appendix  F. 

Other  constimtive  relations  are  given  in  Equations  (3.8-15).  The  dynamic  compaction 
equation  is  expressed  in  Equation  (3.8).  Constituent  one  is  a  gas  described  by  a  virial 
equation  of  state  (3.9).  Constituent  two  is  a  solid  described  by  a  Tait  equation  of  state  [69] 
(3.12).  Assumption  of  a  constant  specific  heat  at  constant  volume  for  each  phase  allows 
caloric  equations  (3.10,13)  and  sound  speed  equations  (3.11,14)  consistent  with  the 
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Figure  3. 1  Comparison  of  Drag  Model  Used  by  Baer  and 
Nunziato  [13]  to  the  Model  of  This  Work 


oil 


r  (pm) 

Figure  3.2  Comparison  of  Heat  Transfer  Model 
Used  by  Baer  and  Nunziato  [13]  to  the 
Model  of  This  Work 
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assumptions  of  classical  thermodynamics  to  be  written  for  each  phase.  Appendix  B  shows 
how  thermodynamically  consistent  equations  are  derived  and  how  relevant  thennodynamic 
properties  are  determined  for  the  state  equations  chosen  here.  The  variable  (()  is  defined  as  a 
volume  fraction,  (j>  =  constituent  volume/total  volume.  Equation  (3.15)  states  that  all  the 
volume  is  occupied  by  constituent  one  or  two;  no  voids  are  permitted. 

By  writing  Equations  (3.1-15)  in  characteristic  form,  it  is  easy  to  show  that  the  model 
is  hyperbolic  and  the  characteristic  wave  speeds  are  u^,  U2^  u^  ±  c^^  and  U2  ±  C2  (see 
Appendix  A).  Baer  has  reached  a  similar  conclusion.  The  fact  that  the  characteristic  wave 
speeds  are  real  is  a  consequence  of  the  assumed  form  of  the  compaction  equation.  Other 
closure  techniques  will,  in  general,  result  in  a  model  with  imaginary  characteristics  which  is 
not  well-posed  for  initial  value  problems. 

The  momentum  and  energy  equations  of  this  model  are  slightly  different  from  those  of 
Baer  and  Nunziato's  model.  The  momentum  equations  of  this  work,  which  are  of  the  same 
general  form  of  those  of  Ref.  1,  differ  with  those  of  Ref.  2  by  a  term  Pj3())j/9x.  Also  the 
energy  equation  of  Ref.  2  includes  a  term  called  "compaction  work,”  proportional  to  the 
volume  fraction  gradient  which  is  not  included  in  this  model.  Which  form  is  correct  is  still 
controversial;  a  defense  of  the  model  presented  here  is  described  in  detail  in  Appendix  C. 
The  methodology  which  is  used  here  to  determine  detonation  structure  is  unaffected  by  the 
particular  choice  of  model  form. 
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IV.  Steady  STATE  COMPACTION  WAVE  Analysis 


This  chapter  is  concerned  with  steady  compaction  waves  in  granular  materials.  These 
waves  are  inert  and  thus  fundamentally  different  from  detonation  waves.  Before  turning  to 
the  study  of  detonation  waves  in  Chapter  5,  there  is  a  good  reason  to  first  consider 
compaction  waves.  That  is,  the  simplicity  of  the  two-phase  equations  allows  a  well- 
understood  solution  to  be  determined.  The  properties  of  this  solution  and  the  solution 
procedure  itself  are  useful  in  the  detonation  analysis. 

A  compaction  wave  can  arise  from  the  impact  of  a  piston  on  a  granular  material.  It  is 
shown  here  that  the  two-phase  equations  are  able  to  describe  such  waves  when  no  reaction 
is  allowed  and  gas  density  is  small  relative  to  the  solid.  This  chapter  has  a  self-contained, 
complete  discussion  of  compaction  waves,  essentially  independent  of  the  detonation 
analysis,  except  that  the  same  model  equations  are  used  in  different  limits.  A  slightly 
different  notation  is  introduced  for  this  chapter  which  reflects  the  simpler  nature  of  the 
compaction  wave  problem  relative  to  the  detonation  wave  problem. 

It  has  been  established  by  experiments  with  granular  high  energy  solid  propellants  [23, 
24]  and  by  numerical  solution  of  unsteady  two-phase  reactive  flow  models  [1,2]  that 
deflagration  to  detonation  transition  (DDT)  in  a  confined  column  of  such  granular  energetic 
material  involves  material  compaction  and  heat  release.  In  many  cases  the  origin  of  such  a 
DDT  can  be  traced  to  the  influence  of  a  compaction  wave,  defined  as  a  propagating 
compressive  disturbance  of  the  solid  volume  fraction  of  the  granular  material.  Steady 
compaction  waves  in  porous  HMX  (cyclic  nitramine)  were  observed  by  Sandusky  and 
Liddiard  [43]  and  Sandusky  and  Bemecker  [44]  arising  from  the  impact  of  a  constant 
velocity  piston  (piston  velocity  <  300  m/s).  Compaction  waves  in  these  experiments  travel 
at  speeds  less  than  800  m/s,  well  below  the  ambient  solid  sound  speed,  which  is  near  3000 
m/s.  To  understand  compaction  waves  it  is  necessary  to  explain  why  this  unusual  result  is 
obtained. 

The  first  step  in  modeling  compaction  waves  is  to  study  steady  compaction  waves. 
With  understanding  gained  from  steady  compaction  waves,  it  is  easier  to  understand  the 
time-dependent  development  of  these  waves  and  how  such  a  wave  can  evolve  into  a 
detonation  wave.  Although  it  is  possible  to  numerically  solve  the  coupled  unsteady  partial 
differential  equations  which  model  such  dynamic  compaction  processes  (including  the 
formation  of  shock  waves)  [56],  it  is  difficult  to  interpret  from  such  models  what  physical 
properties  dictate  the  speed,  pressure  changes,  and  porosity  changes  of  compaction  waves. 
It  is  the  goal  of  this  chapter  to  provide  a  simple  method  to  predict  these  parameters  as  a 
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function  of  material  properties  with  a  representative  model. 

The  experiments  of  Sandusky  and  Liddiard  arc  simulated  by  studying  steady  solutions 
of  two-phase  flow  model  equations.  Without  considering  wave  stmcture,  Kooker  [57]  has 
used  an  algebraic  end  state  analysis  to  predict  compaction  wave  speed  as  a  function  of 
piston  velocity  using  full  two-phase  model  equations.  It  is  possible  to  extend  this  analysis 
in  the  limit  where  the  effect  of  one  of  the  phases  is  dominant  This  approach  was  fu^t  used 
by  Baer  [45]  in  his  study  of  steady  compaction  wave  structure.  Here  a  more  detailed 
discussion  is  provided  of  steady  structure  and  basic  parameter  dependencies.  Throughout 
this  chapter,  the  assumptions  and  results  will  be  compared  to  those  of  Baer. 

The  results  show  a  continuous  dependence  of  compaction  wave  structure  on  the  piston 
velocity  supporting  the  wave;  depending  on  the  piston  velocity,  two  broad  classes  of 
compaction  zone  structures  exist.  At  low  piston  velocities  the  compaction  wave  travels  at 
speeds  less  than  the  ambient  sound  speed  of  the  solid.  Such  waves  are  called  subsonic 
compaction  waves.  The  structure  is  characterized  by  a  smooth  rise  in  pressure  from  the 
ambient  to  a  higher  pressure  equal  to  the  static  pore  collapse  stress  level.  Subsonic 
compaction  waves  have  been  observed  experimentally  (though  compaction  zone  widths 
have  not  been  measured)  and  predicted  by  Baer.  Above  a  critical  piston  velocity  the 
compaction  wave  travels  at  speeds  greater  than  the  ambient  sound  speed  in  the  solid.  A 
discontinuous  shock  wave  leads  a  relaxation  zone  where  the  pressure  adjusts  to  its 
equilibrium  static  pore  collapse  value.  Such  waves  are  called  supersonic  compaction 
waves.  Supersonic  compaction  waves  with  leading  shocks  have  not  as  yet  been  observed 
nor  predicted. 

A  shock  wave  in  compaction  wave  structure  is  admitted  because  the  model  equations 
are  hyperbolic.  This  model  ignores  the  effects  of  diffusive  momentum  and  energy 
transport.  If  included,  these  effects  would  define  the  width  of  the  shock  structure.  Here  it 
is  assumed  that  the  length  scales  on  which  these  processes  are  important  are  much  smaller 
than  the  relaxation  scales  which  define  compaction  zone  structure. 

Compaction  wave  phenomena  predicted  here  have  analogies  in  gas  dynamics.  As 
described  by  Becker  and  Bohme  [58],  gas  dynamic  models  which  include  thermodynamic 
relaxation  effects  predict  a  dispersed  wave  to  result  from  the  motion  of  a  piston  into  a 
cylinder  of  gas.  Steady  solutions  with  and  without  discontinuous  jumps  are  identified. 
These  solutions  have  features  which  are  similar  to  those  predicted  by  the  compaction  wave 
model. 

Here  comments  are  made  on  the  differences  and  similarities  of  the  original  Baer  study 
and  the  present  study.  Baer's  incompressibility  assumption  has  been  relaxed  to  allow  a 
fully  compressible  solid.  A  complete  characterization  of  compaction  wave  structure  as  a 
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function  of  piston  wave  speed  including  an  analysis  of  the  supersonic  case  is  given  here. 
With  this  analysis  many  new  results  are  obtained.  A  unique  equilibrium  condition, 
determined  algebraically,  is  obtained.  As  Baer  does,  it  is  demonstrated  that  the  problem  of 
determining  compaction  wave  structure  can  be  reduced  to  solving  one  ordinary  differential 
equation  for  volume  fraction.  Other  thermodynamic  quantities  (pressure,  density,  etc.)  are 
algebraic  functions  of  volume  fraction.  An  analytic  solution  in  the  strong  shock  limit  is 
given.  A  term  used  by  Baer  called  "compaction  work"  is  not  included  in  this  model.  As 
shown  in  Appendix  C,  this  term  violates  the  principle  of  energy  conservation. 

Unsteady  Model 

The  two-phase  continuum  mixture  equations  (3.1-15)  are  repeated  in  a  condensed 
form  in  Equations  (4.1-7).  The  model  describes  two-phase  flow  with  inter-phase  mass, 
momentum,  and  energy  transport.  A  density,  pi;  pressure.  Pi;  energy,  ei;  temperature,  Tj; 
velocity,  u,;  and  volume  fraction,  is  defined  for  each  phase  (for  the  gas  i  =  1,  for  the 
solid  i  =  2).  A  compaction  equation  similar  to  that  of  Baer  and  Nunziato  is  utilized.  The 
compaction  equation  models  the  time-dependent  pore  collapse  of  a  porous  matrix  and  is 
based  on  the  dynamic  pore  collapse  theory  of  Carroll  and  Holt  [59]. 

The  unsteady  two-phase  equations  are 


KpA)  -  KpA“i)  =  A, 

(4.1) 

KPM)+|^Pi4>i^Pi1>ft)  =  B. 

(4.2) 

|(p,$.[e,  +  ufn])  +  ||p.4,.„[e,  +  uf/2  +  P/p.])  =  c, 

(4.3) 

3<()  0<t)  \ 

(4.4) 

P.  =  P.(p.,  T.) 

1  I  ^1  r 

(4.5) 

e.  =  e.(P.,  p.) 

1  1  1  *^1 

(4.6) 

<j)  +0  =1 
I  2 


(4.7) 
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Equations  (4.1),  (4.2),  and  (4.3)  describe  the  evolution  of  mass,  momentum,  and  energy, 
respectively,  of  each  phase.  Inter-phase  transport  is  represented  in  these  equations  by  the 
terms  Ai,  Bi,  and  Q,  which  are  assumed  to  be  algebraic  functions  of  Pi,  Ui,  pi,  etc.  These 
terms  are  specified  such  that  the  following  conditions  hold: 


(4.8) 


This  insures  that  the  mixture  equations  obtained  by  adding  the  constituent  mass, 
momentum,  and  energy  equations  are  conservative. 

For  each  phase  an  initial  temperature,  density,  velocity,  and  volume  fraction  is 
defined.  The  subscript  0  is  taken  to  represent  an  initial  condition. 

p,  =  p„,  u,„  =  0.  =  (4.9) 

Other  variables  are  determined  by  the  algebraic  relations  (4.5),  (4.6),  and  (4.7). 

Equation  (4.4)  is  the  compaction  equation.  A  sintiilar  model  equation  has  been  used  by 
Butcher,  Carroll,  and  Holt  [60]  to  describe  time-dependent  (dynamic)  pore  collapse  in 
porous  aluminum.  The  parameter  is  defined  as  compaction  viscosity,  not  to  be  confused 
with  the  viscosity  associated  with  momentum  diffusion.  The  compaction  viscosity  defines 
the  only  length  scale  in  this  problem.  The  existence  of  such  a  parameter  is  still  a  modeling 
assumption  and  its  value  has  not  been  determined.  There  is,  however,  a  strong  theoretical 
justification  for  the  dynamic  pore  collapse  model.  It  has  been  shown  (Appendix  A)  that 
when  dynamic  compaction  is  incorporated  into  two-phase  model  equations,  the  equations 
are  hyperbolic.  The  initial  value  problem  is  required  to  be  hyperbolic  in  order  to  insure  a 
stable  solution. 

In  the  compaction  equation  (4.4)  f  represents  the  intra-granular  stress  in  the  porous 
medium.  It  is  assumed  to  be  a  function  of  the  volume  fraction.  Baer  has  estimated  f  from 
Elban  and  Chiarito's  [61]  empirical  quasi-static  data  obtained  by  measuring  the  static 
pressure  necessary  to  compact  a  porous  media  to  a  given  volume  fraction.  Carroll  and  Holt 
have  suggested  an  analytical  form  for  f  for  three  regimes  of  pore  collapse,  an  elastic  phase, 
an  elastic-plastic  phase,  and  a  plastic  phase.  In  this  chapter  f  will  be  modelled  with  an 
equation  similar  to  Carroll  and  Holt’s  plastic  phase  equation.  Here,  two  a  priori 
assumptions  about  f  are  made.  First,  it  is  assumed  that  f  is  a  monotonically  increasing 
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function  of  volume  fraction  so  that  an  increasing  hydrostatic  stress  is  necessary  to  balance 
the  increased  intra-granular  stress  which  arises  due  to  an  increasing  solid  volume  fraction. 
Second,  it  is  assumed  that  at  the  initial  state  f  must  equal  the  difference  of  the  solid  and  gas 
pressures  so  that  the  system  is  initially  in  equilibrium.  The  results  show  that  with  these 
assumptions,  compaction  wave  phenomena  are  relatively  insensitive  to  the  particular 
functional  form  of  f. 

Equations  (4.5)  and  (4.6)  are  state  relations  for  each  phase.  Equation  (4.7)  arises 
from  the  definition  of  volume  fraction.  It  states  that  all  volume  is  occupied  by  either  solid 
or  gas. 

Dimensionless  Steady  Model 

To  study  compaction  waves  in  the  context  of  this  model,  the  following  assumptions 
are  made:  1)  a  steady  wave  travelling  at  speed  D  exists,  2)  gas  phase  equations  may  be 
neglected,  3)  inter-phase  transport  terms  may  be  neglected,  and  4)  the  solid  phase  is 
described  by  a  Tait  equation  of  state.  As  a  result  of  Assumption  1,  Equations  (4.1)  through 
(4.4)  may  be  transformed  to  ordinary  differential  equations  under  the  Galilean 
transformation  ^  =  x  -  Dt,  v  =  u  -  D.  By  examining  the  dimensionless  form  of  Equations 
(4.1)  through  (4.7),  it  can  be  shown  that  in  the  limit  as  the  ratio  of  initial  gas  density  to 
initial  solid  density  goes  to  zero,  that  there  is  justification  in  neglecting  gas  phase  equations 
and  inter-phase  transport.  To  prove  this  contention,  one  can  integrate  the  steady  mixture 
mass,  momentum,  and  energy  equations  formed  by  adding  the  component  equations  to 
form  algebraic  mixture  equations.  By  making  these  equations  dimensionless  (as  done  in 
Chapter  5),  it  is  seen  that  all  gas  phase  quantities  are  multiplied  by  the  density  ratio 
Pio/p2o-  long  as  dimensionless  gas  phase  properties  are  less  than  O(p2o/Pio)>  there  is 
justification  in  neglecting  the  effect  of  the  gas  phase. 

Because  the  gas  phase  is  neglected,  the  subscripts  1  and  2  are  discarded.  All  variables 
are  understood  to  represent  solid  phase  variables.  The  caloric  Tait  equation  [69]  for  the 
solid  is 


P  +  p  s 

e  - - a-  (4.10) 

(Y-  Dp 


Here  y  and  s  are  parameters  that  define  the  Tait  state  equation.  The  value  of  y  is  chosen  to 
match  shock  Hugoniot  data  [55].  It  is  analogous  to  the  specific  heat  ratio  for  an  ideal 
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equation  of  state.  The  parameter  s  is  defined  as  the  non-ideal  solid  parameter.  In  this  study 
s  is  viewed  as  an  adjustable  parameter  which  allows  the  equation  of  state  to  be  varied  in  a 
simple  way  in  order  to  show  how  the  results  are  sensitive  to  non-ideal  state  effects.  When 
s  =  0,  the  state  equation  is  an  ideal  state  equation.  For  this  study  a  value  of  s  was  chosen  to 
match  the  compaction  wave  data  of  Sandusky  and  Uddiard  [43]. 

To  determine  the  ambient  solid  sound  speed,  an  important  term  in  this  analysis,  it  is 
necessary  to  specify  a  thermal  equation  of  state.  By  assuming  a  constant  specific  heat  at 
constant  volume  c^,  a  thermal  equation  of  state  consistent  with  Equation  (4.10)  can  be 
derived. 


P  =  (v-Dc^pT  -  p^s/Y  (4.11) 

Based  on  Equations  (4.10)  and  (4.1 1)  an  equation  for  the  solid  sound  speed  c  is  easily 
derived  by  using  the  thermodynamic  identity  T  dt)  =  de  -  P/p2  dp,  where  tj  is  the  entropy. 

=  ~L  =  Y(Y-1)c^T  (4.12) 

op  ' 

To  simplify  the  analysis,  dimensionless  variables  are  denoted  by  a  star  subscript  and 
are  defined  as  follows 


P*  ^  ^  ^0  ’  ~  ^  f  ^  f  e^  —  e  /  D  ,  T*  =  c^T  /  D  , 

P.  =  p/(py),  ^  =  5p^d/p^ 


With  this  choice  of  dimensionless  variables  four  dimensionless  parameters  arise. 


Y  =  Tait  Solid  Parameter,  - 

2 

D  Y 


=  a  =  Non-Ideal  Solid  Parameter 


Y(Y-  ^ 


YD^ 


=  7C  =  initial  pressure  ,  (J)^  =  initial  volume  fraction 


25 


For  materials  of  interest  <l)o  and  y  are  of  order  1.  Interesting  limiting  cases  can  be  studied 
when  s— >  0,  corresponding  to  either  the  strong  shock  or  weak  non-ideal  effect  limit,  or 
when  7C  — >  0,  corresponding  to  the  strong  shock  limit 

With  the  assumptions  made,  steady  dimensionless  equations  can  be  written  to  describe 
the  compaction  of  an  inert  solid  porous  material  as  follows: 


dq* 

(4.13) 

— (p*<!>  +  P.<}>v*)  =0 

(4.14) 

-( 

P*<l>v*^e*  + vJ/2  +  P*/p,j  j  =  0 

(4.15) 

(4.16) 

P^+ya 

(4.17) 

(Y-1)  P* 


Initial  conditions  are  specified  as 

=  1  ,  <1>  =  <!>q  ,  =  -1  ,  =  7t  (4.18) 

Equations  (4.13-17)  are  equivalent  to  Baer's  steady  model  except  a  term  Baer  calls 
"compaction  work"  is  not  included  and  a  simpler  state  equation  is  used.  Equations  (4.13), 
(4.14),  and  (4.15)  may  be  integrated  subject  to  initial  conditions  (4.18)  resulting  in  the 
following  set  of  equations: 


P*  -  f,  ((|>)) 


P*  0  V*  =  -  A 


(4.20) 
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P*<l>  +  P*<l>vJ  =  <{)J 

ir  +  1) 

(4.21) 

1  n  +  yo  \ 

P*<l)vJe*  +  v*/2  +  P^/p*l  =  .(|) 

- l/2  +  7t 

(4.22) 

\  /  0 

\  Y-l  j 

P*+YC 

e*  - - 

(4.23) 

(Y- 1)  P* 

From  Equations  (4.20)  through  (4.23),  equations  for  pressure  and  velocity  as 
functions  of  volume  fraction  can  be  written.  Equation  (4.23)  is  used  to  eliminate  energy 
from  Equation  (4.22).  Velocity  is  eliminated  from  Equations  (4.21)  and  (4.22)  by  using 
Equation  (4.20).  Then  density  is  eliminated  from  Equation  (4.22)  by  using  Equation 
(4.21).  What  remains  is  a  quadratic  equation  involving  only  pressure  and  volume  fraction. 
It  is  possible  to  solve  this  equation  for  pressure  explicitly  in  terms  of  volume  fraction.  The 
solution  is 


P* 


"’o 

■KV+dI  '*’0  I 


-1  ± 


(4.24) 


The  solution  corresponding  to  the  positive  branch  is  the  physically  relevant  one.  The 
negative  branch  is  associated  with  negative  pressure.  Equations  (4.20)  and  (4.21)  may  be 
simultaneously  solved  for  velocity  as  a  function  of  pressure  and  volume  fraction.  The 
velocity  is  given  by 


V*  = 


P^  »  -  (1  +  7t) 


(4.25) 
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By  using  Equation  (4.24)  to  substitute  for  pressure  in  Equation  (4.25),  velocity  is  available 
as  a  function  of  volume  fraction  alone.  The  mass  equation  (4.20)  can  be  used  to  give 
density  as  a  function  of  volume  fraction  and  then  the  state  equation  (4.17)  can  be  used  to 
give  energy  as  a  function  of  volume  fraction.  Thus  all  variables  in  the  compaction  equation 
(4.19)  can  be  expressed  as  functions  of  volume  fraction;  the  compaction  wave  problem  is 
reduced  to  solving  one  ordinary  differential  equation  (4.19)  for  volume  fraction  subject  to 
the  condition  ([)  =  (t)o  at  =  0. 

Next  the  technique  is  described  for  determining  wave  speed  as  a  function  of  piston 
velocity.  This  calculation  is  algebraic  and  can  be  made  without  regards  to  structure.  The 
solution  is  parameterized  by  the  wave  velocity  through  the  definitions  of  k  and  a.  Instead 
of  using  a  piston  velocity  as  an  input  condition,  it  is  easier  to  consider  the  wave  speed  to  be 
known  and  from  that  wave  speed  calculate  a  piston  velocity.  By  assuming  a  static  pressure 
equilibrium  end  state  in  Equation  (4.19)  (?•(<())  =  f*((t>)),  it  is  possible  to  determine  the 
equilibrium  volume  fraction  and  thus,  from  Equations  (4.24)  and  (4.25),  the  final  velocity 
V*.  The  piston  velocity  (Up)  is  found  by  transforming  the  final  velocity  to  the  lab  frame  by 
using  the  transformation  Up  =  D(v*  +  1). 

Pressure  equilibrium  end  states  are  found  when  a  volume  fraction  is  found  such  that 
the  pressure  given  by  Equation  (4.24)  matches  the  intra-granular  stress  predicted  by  f.  In 
the  initial  state.  Equation  (4.24)  predicts  a  pressure  of  re,  the  dimensionless  initial  pressure. 
By  assumption  f  also  yields  a  value  of  k  in  the  initial  state  so  that  the  undisturbed  material 
is  stationary.  In  Figure  4.1,  dimensional  pressure  in  HMX  is  plotted  as  a  function  of 
volume  fraction  from  Equation  (4.24)  for  a  series  of  wave  speeds  and  an  initial  volume 
fraction  of  0.73.  Except  for  compaction  viscosity  Pp  parameters  used  to  model  HMX  are 
those  previously  listed  in  Table  I  (c^,  lip,  and  Po  of  Baer  is  used,  and  the  parameters  y  and 
s  are  estimated  by  requiring  predictions  to  match  shock  and  compaction  data.  Unlike  in 
detonation  wave  analysis,  there  is  no  special  problem  posed  by  using  Baer’s  value  of 
compaction  viscosity,  1000  kg/(m  s),  in  these  calculations).  All  curves  pass  through  the 
point  of  initial  pressure  and  volume  fraction. 

The  curve  on  Figure  4.1  for  the  ambient  sonic  wave  speed  (D  =  3000  m/s)  has  a 
special  property  whose  importance  will  be  apparent  in  the  following  discussion.  For  this 
curve  a  volume  fraction  minimum  exists  at  the  initial  volume  value.  It  can  be  proven  for  a 
sonic  wave  speed,  that  the  discriminant  in  Equation  (4.24)  is  identically  zero  for  <J)  =  <t>o  and 
D  =  y(y  -  Uc^Tq  (the  ambient  solid  sonic  wave  speed). 

The  positive  pressure  branch  of  Equation  (4.24)  is  a  double- valued  function  of  volume 
fr-action  for  wave  velocities  that  exceed  the  ambient  solid  sound  speed  and  single-valued  for 
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Figure  4.1  Pressure  vs.  Volume  Fraction  for  Subsonic,  Sonic,  and  Supersonic  Wave  Speeds 


wave  velocities  less  than  or  equal  to  the  ambient  solid  sound  speed.  For  subsonic  wave 
speeds,  small  increases  from  the  initial  volume  fraction  cause  small  positive  perturbations 
in  pressure.  For  supersonic  wave  speeds  a  positive  increase  of  the  initial  volume  fraction  is 
only  acceptable  if  the  pressure  jumps  discontinuously  to  a  shocked  value  on  the  upper 
portion  of  the  double-valued  P*-<|)  curve.  Because  the  governing  equations  are  hyperbolic, 
these  shock  jumps  are  admissible.  From  Equation  (4.19)  the  shock  jump  condition  for 
volume  fraction  is 


Ul-o 


(4.26) 


where  "0"  denotes  the  initial  state  and  "s"  the  shock  state.  Thus  the  shock  volume  fraction 
is  always  equal  to  the  initial  volume  fraction. 

From  Equations  (4.24)  and  (4.25)  the  shock  pressure  and  particle  velocity  can  be 
determined.  The  shocked  values  are  independent  of  the  initial  solid  volume  fraction. 


P  = 


2-(jc  +  a)(Y- 1) 


-  a 


V  =  - 


Y+  1 

(Y-  l)  +  2Y(Jt  +  <y) 
Y  +  1 


(4.27) 

(4.28) 


The  combination  of  parameters  tc  -h  a  is  independent  of  the  non-ideal  solid  parameter  s.  So 
from  Equations  (4.27)  and  (4.28)  it  is  deduced  that  non-ideal  effects  lower  the  shock 
pressure  by  a  constant,  o,  and  do  not  affect  the  shock  particle  velocity. 

Based  on  the  implications  of  Equation  (4.24),  the  structure  analysis  is  thus 
conveniently  split  in  two  classes,  subsonic  and  supersonic.  As  wave  speed  increases  from 
subsonic  values,  the  initial  pressure  at  the  wave  front  is  the  ambient  pressure  until  the 
compaction  wave  speed  is  sonic.  For  wave  speeds  greater  than  the  ambient  solid  sound 
speed  the  initial  pressure  jumps  are  dictated  by  Equation  (4.27).  A  plot  of  the  leading 
pressure  versus  compaction  wave  speed  is  shown  in  Figure  4.2. 
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Figure  4.2  Pressure  at  Compaction  Wave  Head  vs.  Compaction  Wave  Speed 

As  an  aside,  it  is  noted  that  a  criterion  fOT  a  solid  equation  of  state  is  that  the  candidate 
equation  along  with  the  Rankine-Hugoniot  jump  conditions  be  able  to  match  experimental 
piston  impact  data.  Typically  parameters  for  solid  equations  of  state  are  determined  by 
choosing  them  such  that  shock  data  is  matched.  For  voidless  HMX  (<j)  =  1)  observations  of 
shock  wave  speed  as  a  function  of  piston  velocity  are  reported  by  Marsh  [55].  By 
rewriting  Equation  (4.28)  in  dimensional  form,  the  wave  speed  D  is  solved  for  as  a 
function  of  piston  velocity. 


D  = 


+ 


(4.29) 


From  Equation  (4.12),  the  term  Y  (y-l)  Cy  To  is  the  square  of  the  ambient  sound  speed  for 
the  non-ideal  solid.  In  a  result  familiar  from  gas  dynamics,  it  can  be  deduced  from 
Equation  (4.29)  that  the  minimum  steady  shock  wave  speed  admitted  in  response  to  a 
piston  boundary  condition  is  the  ambient  sonic  speed.  For  values  of  y,  Cy,  and  Tq  listed  in 
Table  I,  the  shock  wave  speed  D  is  plotted  as  a  function  of  piston  velocity  Up  and  data  from 
Marsh  in  Figure  4.3. 
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Figure  4.3  Piston  Velocity  vs.  Solid  Shock  Speed 

The  parameter  y  has  been  fixed  such  that  there  is  agreement  between  the  data  and  the  model 
predictions.  In  the  range  of  piston  velocities  shown.  Equation  (4.29)  approximates  a  linear 
D  vs.  Up  relation  used  by  other  modelers  to  match  this  data. 

Subsonic  Compaction  Waves 


To  study  subsonic  compaction  waves  admitted  by  Equation  (4.19),  a  form  for  f*  is 
chosen: 


This  function  satisfies  the  requirements  described  earlier,  namely,  it  is  a  monotonically 
increasing  function  of  volume  fraction  and  is  constructed  such  that  the  system  is  in 
equilibrium  in  the  initial  state.  It  has  the  same  form  as  the  plastic-phase  static  pore  collapse 
relation  given  by  Carroll  and  Holt  [59].  It  is  not  the  Carroll  and  Holt  relation,  as  the 
leading  coefficient  in  the  Carroll  and  Holt  relation  is  the  yield  stress  of  the  solid.  In 


Equation  (4.30)  the  leading  coefficient  is  a  function  of  initial  volume  fraction.  Predictions 
of  Equation  (4.30)  approximately  match  the  experimental  results  of  Elban  and  Chiarito 
[61].  Figure  4.4  compares  a  curve  fit  of  Elban  and  Chiarito's  data  for  HMX  with  the 
approximation  given  by  Equation  (4.30). 


Equation  (4.30) 


Figure  4.4  Comparison  of  Static  Pore  Collapse  Data  with  Predictions  of  Equation  (4.30) 

To  locate  an  end  state.  Equations  (4.24)  and  (4.30)  are  solved  simultaneously.  For 
73%  theoretical  maximum  density  (TMD)  HMX  (volume  fraction  =  0.73)  and  a  variety  of 
subsonic  wave  speeds,  curves  of  pressure  versus  volume  fraction  from  Equations  (4.24) 
and  (4.30)  are  plotted  in  Figure  4.5.  As  wave  speed  increases,  the  final  volume  fraction 
increases.  For  wave  speeds  above  600  m/s  nearly  complete  compaction  is  predicted.  For 
wave  speeds  of  about  200  m/s  or  lower,  no  steady  compaction  wave  is  predicted.  This  is 
solely  a  consequence  of  the  assumed  form  of  f.  The  form  of  f  chosen  crosses  through  the 
initial  point  with  a  positive  slope  and  fails  to  intersect  the  pressure-volume  fraction  curves 
for  low  wave  speeds. 

For  73%  TMD  HMX  Figure  4.6  shows  plots  of  compaction  wave  speed,  final 
density,  final  volume  fraction,  final  pressure,  and  final  mixture  pressure  (mixture  pressure 
=  pressure  *  volume  fraction)  versus  piston  velocity.  Also  shown  are  the  observations  of 
Sandusky  and  Liddiard  [43]  and  Sandusky  and  Bemecker  [44]  of  wave  speed  and  final 
volume  fraction  and  their  predictions  of  pressure.  The  relatively  small  density  changes 
verify  that  Baer's  incompressibility  assumption  is  a  good  approximation.  Figure  4.7 
shows  predictions  of  compaction  wave  speed,  final  volume  fraction,  and  final  mixture 
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Figure  4.5  Pressure  vs.  Volume  Fraction  for  Subsonic  Wave  Speeds  and  Pore  Collapse  Function,  f 
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Figure  4.6  Compaction  Wave  End  States  vs.  Piston  Velocity 


36 


pressure  as  a  function  of  initial  volume  fraction  for  a  constant  piston  velocity  of  100  m/s 
along  with  Sandusky's  predictions  as  reported  by  Kooker  [62]. 

Subsonic  Structure 


Equation  (4.19)  has  been  numerically  integrated  to  determine  the  structure  of  the 
subsonic  compaction  zone.  The  integration  was  performed  using  the  IMSL  routine 
DVERK,  a  fifth  and  sixth  order  Runge-Kutta  routine.  A  step  size  was  chosen  such  that  the 
compaction  zone  stmcture  was  described  by  about  100  points.  Using  more  points  had  little 
effect  on  the  results.  Run  times  to  determine  a  structure  were  less  than  ten  seconds  on  the 
UIUC  Cyber  175  computer.  In  the  numerical  integrations  pressure,  velocity,  and  f  are 
used  as  given  by  Equations  (4.24),  (4.25),  and  (4.30),  respectively.  The  integration  was 
performed  starting  at  ^*  =  0  and  integrating  towards  -4  -o®.  To  initiate  the  integration,  a 
small  positive  perturbation  of  volume  fraction  was  introduced  which  in  this  case  causes  a 
small  positive  perturbation  in  pressure. 

Figure  4.8  shows  the  particle  velocity,  volume  fraction,  and  pressure  in  the 
compaction  zone  for  a  subsonic  compaction  wave.  Here  the  piston  velocity  is  100  m/s  and 
the  initial  volume  fraction  is  0.73.  The  compaction  wave  speed  is  404.63  m/s.  For  an 
assumed  compaction  viscosity  of  1000  kg/(m  s)  a  compaction  wave  thickness  of  80  mm  is 
predicted.  Because  compaction  viscosity  defines  the  only  length  scale  in  this  problem, 
compaction  viscosity  only  serves  to  define  the  compaction  wave  thickness.  For  the  same 
value  of  compaction  viscosity  Baer  reports  a  compaction  wave  thickness  of  31.9  mm.  The 
discrepancy  could  be  due  to  many  effects  including  the  definition  of  compaction  zone 
length.  It  is  important  to  note  that  the  length  is  of  the  same  order  of  magnitude.  Final 
pressure,  wave  speed,  and  final  volume  fraction  are  unaffected  by  the  value  chosen  for 
compaction  viscosity.  By  measuring  a  compaction  wave  thickness,  an  estimate  could  be 
made  for  the  compaction  viscosity. 

Supersonic  Compaction  Waves 

Supersonic  End  States 

At  0.73  initial  porosity  for  piston  velocities  greater  than  884  m/s,  supersonic 
compaction  waves  are  also  admitted.  Figure  4.9  shows  plots  of  compaction  wave  speed, 
final  density,  final  volume  fraction,  final  pressure,  and  final  mixture  pressure  as  a  function 
of  piston  velocity.  These  curves  encompass  both  the  subsonic  and  supersonic  compaction 
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Figure  4.9  Subsonic  and  Supersonic  Compaction  End  States  vs.  Piston  Velocity 
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wave  end  states.  It  is  seen  that  the  end  states  are  a  continuous  function  of  piston  velocity. 
In  Figure  4.9  the  shock  wave  speed  as  a  function  of  piston  velocity  is  plotted  alongside  the 
compaction  wave  speed.  For  large  wave  speeds  the  predicted  shock  velocity  converges 
with  the  compaction  wave  velocity.  It  is  demonstrated  next  that  this  is  a  consequence  of 
non-ideal  effects  having  litde  importance  at  supersonic  wave  speeds.  Furthermore  it  will  be 
demonstrated  that  the  existence  of  subsonic  compaction  waves  can  be  attributed  solely  to 
non-ideal  effects. 

Supersonic  Structure 

Equations  (4.24)  and  (4.25)  can  be  simplified  in  the  limit  as  a  0.  The  limit  of  small 
a  corresponds  either  to  negligible  non-ideal  effects  or  large  wave  speed.  In  the  limit  as  a 
0  Equations  (4.24),  (4.25),  and  (4.19)  can  be  written  as 


P* 


V 


V 

s 


(4.31) 

(4.32) 

(4.33) 


Equation  (4.32)  holds  that  in  this  limit  the  velocity  is  constant  in  the  relaxation  zone 
and  is  equal  to  the  shocked  particle  velocity.  For  s  s  0  (that  is  for  an  ideal  state  relation) 
Equation  (4.32)  is  equivalent  to  Equation  (4.29);  thus,  for  an  ideal  state  relation  the 
minimum  compaction  wave  speed  is  the  ambient  sonic  speed.  Any  subsonic  compaction 
wave  admitted  by  the  model  (Equations  (4.19)  -  (4.23))  is  a  direct  consequence  of  non¬ 
ideal  state  effects. 

In  the  strong  shock  limit  D  -» <»,  and  both  k  and  ct  0.  Equation  (4.33)  has  a 
simple  solution  in  this  limit,  assuming  f*  to  be  sufficiently  bounded.  (Note  that  because  of 
the  logarithmic  singularity  at  (()  =  1,  that  Equation  (4.30)  does  not  meet  this  criterion.  The 
general  model  is  not,  however,  restricted  to  a  function  of  this  form)  In  this  limit  Equation 
(4.33)  becomes 
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d(j> 


whose  solution  is 


(4.34) 


(4.35) 


In  terms  of  dimensional  parameters,  the  compaction  zone  thickness  found  by  equating  the 
exponent  in  Equation  (4.35)  to  one  and  substituting  the  expression  for  piston  velocity  for 
wave  speed  is  estimated  as 


comp 

The  length  is  proportional  to  compaction  viscosity  and  inversely  proponional  to  piston 
velocity  and  the  product  of  density  and  volume  fraction. 

An  example  of  supersonic  structure  arising  from  the  impact  of  a  1000  m/s  piston  is 
now  given.  Figure  4.10  shows  the  particle  velocity,  volume  fraction,  and  pressure  in  the 
compaction  zone  for  a  supersonic  compactiqn  wave.  Here  the  initial  volume  fraction  is 
0.73.  The  compaction  wave  speed  is  3353.67  m/s  and  the  wave  thickness  is  2.9  mm.  It  is 
seen  that  pressure  and  particle  velocity  undergo  shock  jumps.  Volume  fraction  does  not 
undergo  a  shock  jump;  however,  its  derivative  does  jump  at  the  initial  point. 

Compaction  Zone  Thickness 

It  is  possible  to  study  the  parametric  dependence  of  compaction  zone  thickness.  Given 
a  constant  compaction  viscosity,  the  model  can  predict  compaction  zone  thickness  as  a 
function  of  initial  volume  fraction  and  piston  velocity.  Should  experiments  be  devised  to 
measure  the  compaction  zone  thickness,  the  experiments  could  provide  a  means  to  verify 
the  theory. 

The  thickness  is  defined  as  the  distance  at  which  the  ratio  of  the  difference  of 
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Figure  4.10  Supersonic  Compaction  Zone  Structure 
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instantaneous  volume  fraction  and  initial  volume  fraction  to  the  difference  of  final  volume 
fraction  and  initial  volume  fraction  is  equal  to  0.99.  The  final  volume  fraction  is  available 
from  the  algebraic  end  state  calculation.  Figure  4.11a  shows  the  compaction  zone  length 
versus  initial  volume  fraction  for  a  piston  velocity  of  100  m/s  and  compaction  viscosity  of 
1000  kg/(m  s).  It  is  not  understood  why  the  peak  in  this  curve  occurs.  It  is  noted  that  for 
high  initial  volume  fraction,  the  zone  length  decreases  as  initial  volume  fraction  increases  in 
accordance  with  the  predictions  of  Equation  (4.36)  for  supersonic  compaction.  It  is 
speculated  that  for  low  porosity  a  different  mechanism  dictates  the  subsonic  compaction 
zone  length  t^’an  supersonic  length.  Figure  4.11b  shows  compaction  zone  length  versus 
piston  velocity  for  73%  TMD  HMX  and  compaction  viscosity  of  1000  kg/(m  s).  The 
compaction  zone  length  decreases  with  increasing  piston  velocity  in  accordance  with  the 
predictions  of  Equation  (4.36)  for  supersonic  compaction.  Figure  4.11c  shows 
compaction  zone  length  as  a  function  of  compaction  viscosity  for  a  100  m/s  piston  velocity 
and  0.73  initial  volume  fraction.  As  no  estimates  are  available  for  compaction  viscosity, 
compaction  zone  lengths  for  a  wide  range  of  compaction  viscosity  have  been  plotted. 
Though  plotted  on  a  log  scale,  the  relationship  is  truly  linear  with  the  compaction  zone 
length  equal  to  a  constant  multiplied  by  the  compaction  viscosity. 


Compaction  Zone  Length  (m)  Compaction  Zone  Length  (m)  Compaction  Zone  Length  (m) 


ig.  4.1  la  Compaction  zone  length  vs.  initial  volume  fraction 


Fig.  4. 1  Ic  Compaction  zone  length  vs.  compaction  viscosity 
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V.  Steady  State  Detonation  Wave  analysis 


Equations  (3.1-15)  can  be  re-cast  in  a  more  tractable  form  using  the  steady  state 
assumption.  First  the  equations  are  written  in  dimensionless  form.  For  a  right-running 
steady  wave  the  Galilean  transformation  ^  =  x-Dt,  v  =  u-  D  causes  Equations  (3.1-8)  to 
become  eight  ordinary  differential  equations.  Here  D  is  a  constant  defined  as  the  steady 
wave  speed.  Next,  Equations  (3.1),  (3.3),  and  (3.5)  may  be  eliminated  in  favor  of 
homogeneous  mixture  equations  formed  by  the  addition  of  the  steady  form  of  Equations 
(3.1)  and  (3.2),  (3.3)  and  (3.4),  and  (3.5)  and  (3.6),  respectively.  The  resultant  mixture 
equations  and  the  steady  form  of  Equation  (3.7)  may  be  integrated  to  form  algebraic 
equations.  Thus  the  steady  two-phase  model  is  described  by  four  ordinary  differential 
equations  and  eleven  algebraic  equations. 

Dimensionless  Steady  Equations 

To  reduce  the  number  of  independent  parameters,  dimensionless  equations  are 
introduced.  Define  dimensionless  variables  where  indicates  a  dimensionless  quantity; 

5.=  VL  V.,  =  v./D  P..  =  P./(p,„D=') 

P.i  =  Pi'PiO  S  =  T..=c^T/D^  r.  =  r/L 

i  =  1,2 

Here  D  is  the  wave  speed;  L,  a  length  scale  which  can  be  associated  with  the  reaction  zone 
length;  p.^^,  the  initial  density  of  phase  i;  c^.,  the  specific  heat  at  constant  volume  of  phase  i. 
Define  the  following  independent  dimensionless  parameters  as 

"l  =  ‘  '  "2  =  I*  '  V  (PjO^lD) 

’'4  =  “  ^  =  P^o'P20  %  =  = 

=  p^DL/p^  =  "n'l’io 
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’'12  =  ’'13  '  '’P 


%  =C-T„/D^ 

10  14  v2  0 


=  Y 
17  *2 

and  the  following  dependent  dimensionless  parameters  as 
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Then  the  dimensionless  model  differential  equations  (for  compact  notation  the  stars  are 
dropped)  can  be  written  as 


(5.1) 

(5.2) 

(5.3) 

(5.4) 
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The  supplemental  algebraic  mixture,  number  conservation,  state,  and  volume  fraction 
equations  are 


Pj1>,Vi  +  -P,1>2V2  =  -Tt 
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(5.5) 

(5.6) 

(5.7) 

(5.8) 

(5.9) 

(5.10) 

(5.11) 

(5.12) 

(5.13) 

(5.14) 

(5.15) 


Undisturbed  conditions  for  the  ambient  mixture  are 


p  =  1  V-  =  -  1  0^  =  It, 


20 


^2  ^14 
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Equilibrium  End  State  Analysis 


To  place  a  first  restriction  on  the  steady  solutions  admitted  by  Equations  (5.1-15), 
equilibrium  end  states  are  considered.  It  is  later  shown  that  the  complete  reaction  state  is  an 
equilibrium  state  for  Equations  (5.1-4).  This  result  can  be  used  at  this  point  to  completely 
describe  the  gas  phase  equilibrium  state.  In  the  complete  reaction  state  the  mixture 
equations  (5.5-7)  allow  for  the  gas  phase  properties  to  be  determined.  For  (t)^  =  0  ( (t)^  =  1), 
Equations  (5.5-7)  can  be  combined  to  form  an  equivalent  two-phase  Rayleigh  line  (5.16) 
and  two-phase  Hugoniot  (5.17) 


^llVl4 


7C 


18 


(5.17) 


From  the  state  relations  (5.9,10)  the  energy  ej  can  be  written  as  ei(Pi,Pi)  which  is 
substituted  into  the  Hugoniot  equation  (5.17).  The  Rayleigh  line  equation  (5.16)  allows  pj 
to  be  eliminated  in  favor  of  Pi-  Substituting  this  in  the  reduced  Hugoniot  equation  results 
in  a  cubic  equation  for  P^  Depending  on  the  wave  speed  three  cases  are  possible:  three 
distinct  real  solutions,  two  equal  real  solutions  and  a  third  real  solution,  and  a  real  solution 
and  a  pair  of  complex  conjugate  solutions.  When  three  distinct  real  solutions  exist,  two  are 
analogous  to  the  weak  and  strong  solutions  predicted  by  the  simple  one-phase  theory.  The 
third  solution  has  no  such  counterpart  and  often  is  a  nonphysical  solution  with  Pi  <  0.  A 
sketch  of  the  two-phase  Rayleigh  lines  and  Hugoniots  for  wave  speeds  corresponding  to 
the  three  classes  of  solutions  is  shown  in  Figure  5.1. 


Figure  5. 1  Sketch  of  Two-Phase  Complete  Reaction  Rayleigh  Line  and  Hugoniot 
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By  imposing  the  condition  that  two  real  roots  are  degenerate  (which  forces  the 
Rayleigh  line  and  Hugoniot  to  be  tangent)  a  minimum  detonation  velocity  can  be  found. 
This  will  be  called  the  CJ  condition.  Because  the  detonation  velocity  D  is  contained  in  the 
dimensionless  parameters,  it  is  convenient  to  return  to  dimensional  variables  to  express  CJ 
conditions.  Define  the  bulk  density,  bulk  pressure,  and  bulk  internal  energy  p^,  P^,  and  e^: 


p  =  p  <()  +  p  d) 

^10^10  ^20^20 


P  =  P  d)  +  P  d) 
a  20'''20 


e 

a 


^10^10^10  ^  ^20^20^20 

p  d)  +  p  d) 
^10^10  ^20^20 


(5.18,  19) 
(5.20) 


The  dimensional  equations  which  must  be  solved  to  determine  the  two-phase  CJ  end  state 
are  shown  next. 


Pl-  Pa-^Pa'D'(l/p,-l/p,) 

(P^^P,)(l/p,-l/p^)  ^ 


C  iPi 
vl  1 


RPj  (bpj  +  1) 


-e  =  0 


dp. 


5.21 


dp. 


5.22 


(5.21) 

(5.22) 

(5.23) 


Equation  (5.21)  is  the  dimensional  form  of  the  Rayleigh  line  equation  (5.16),  Equation 
(5.22)  is  the  dimensional  form  of  the  Hugoniot  equation  (5.17),  and  Equation  (5.23)  is  the 
tangency  condition.  These  three  equations  have  been  solved  by  iteration  for  the  three 
unknowns,  Pj,  p^,  and  D.  The  equations  have  an  exact  solution  in  terms  of  a  quadratic 
equation  in  the  ideal  gas  limit  (b  =  0).  The  ideal  gas  solution  has  been  used  as  a  first 
estimate  for  the  iterative  solution. 

The  effects  of  a  non-ideal  gas  and  small  initial  bulk  pressure  on  CJ  conditions  can  be 
seen  by  writing  the  CJ  conditions  as  Taylor  series  expansions  which  are  valid  in  the  limit  as 
the  dimensionless  groups  bp^  and  PaAPa^a)  approach  zero.  These  expressions  were 
obtained  with  the  aid  of  the  computer  algebra  program  MACSYMA  and  have  been  verified 
by  comparing  predictions  to  the  solutions  obtained  by  iteration.  (The  same  technique  can 
be  used  to  obtain  CJ  deflagration  conditions  for  a  two-phase  material;  these  conditions  are 


49 


reported  in  Appendix  D).  The  Taylor  series  expansions  for  the  CJ  detonation  condition  for 
a  two-phase  material  are  given  below.  The  leading  coefficients  on  the  right  hand  sides  of 
Equations  (5.24-28)  are  the  exact  solutions  in  the  limit  of  no  non-ideal  gas  effects  (b  =  0) 
and  zero  initial  bulk  pressure  (P^  =  0).  The  bracketed  terms  in  these  equations  represent  the 
first  order  corrections  for  finite  non-ideal  effects  and  finite  initial  bulk  pressure.  From  the 
expressions,  it  is  seen  that  non-ideal  effects  tend  to  raise  the  detonation  wave  speed  and 
pressure,  lower  the  density,  and  have  no  effect  on  the  gas  velocity  or  temperature.  Finite 
initial  bulk  pressure  tends  to  lower  the  detonation  wave  speed,  pressure,  density,  and  gas 
velocity,  and  has  no  effect  on  the  temperature  at  this  order  of  the  expansion. 
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(5.24) 


(5.25) 


(5.26) 


(5.27) 


(5.28) 


For  =  b  =  0  these  formulae  show  that  it  is  appropriate  to  treat  the  two-phase  CJ 
condition  as  a  one-phase  CJ  condition  using  Pg  and  eg  as  effective  one-phase  properties. 
Fickett  and  Davis  [40]  give  equations  for  one-phase  CJ  properties  for  an  ideal  gas  in  the 
limit  of  small  initial  pressure.  In  these  equations,  one  can  simply  substitute  the  bulk 
density  for  the  initial  density  and  the  bulk  internal  energy  for  the  chemical  energy  to  obtain 
the  two-phase  CJ  equations.  It  is  important  to  note  that  the  two-phase  CJ  properties  are 
predicted  from  the  full  model  equations.  The  two-phase  nature  of  the  conditions  is 
embodied  in  the  definitions  of  bulk  properties,  which  have  no  one-phase  counterpart. 

Figure  5.2  shows  plots  of  the  CJ  properties  predicted  by  Equations  (5.24-28)  along 
with  the  exact  CJ  properties  predicted  by  iterative  solution  of  Equations  (5.21-23).  Also 
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plotted  on  these  curves  are  predictions  of  CJ  properties  from  the  thermochemistry  code 
TIGER  [54]  as  reported  in  Refs.  1  and  47.  It  is  seen  that  the  predictions  of  the 
approximate  formulae  more  accurately  predict  the  exact  solutions  for  low  initial  bulk 
density.  The  improved  accuracy  for  low  initial  bulk  density  can  be  attributed  to  the  form  of 
the  Taylor  series  prediction,  whose  accuracy  improves  as  the  dimensionless  parameter  bp^ 
approaches  zero.  Except  for  the  CJ  density,  the  approximate  formulae  estimate  the  general 
trends  for  a  large  range  of  initial  bulk  densities. 

Equations  (5.24,  25)  indicate  that  the  CJ  state  is  quite  sensitive  to  the  non-ideal 
parameter  b,  a  parameter  allowed  to  vary  in  Ref.  47  to  match  CJ  TIGER  predictions.  In 
particular,  when  the  dimensionless  group  bp^  is  of  order  1,  noh-ideal  effects  become  quite 
important.  This  is  demonstrated  in  Figure  5.3  which  for  constant  bulk  density  plots  CJ 
wave  speed  versus  the  non-ideal  parameter  b.  This  plot  was  obtained  by  solving  the  full 
non-linear  equations  (5.21-23). 


b  (m^  /kg) 


Figure  5.3  CJ  Wave  Speed  vs.  Non-Ideal  Parameter  b 

By  numerically  studying  exact  two-phase  CJ  conditions,  it  can  be  inferred  that  the  CJ 
point  is  a  sonic  point;  that  is,  the  gas  velocity  relative  to  the  wave  head  is  equal  to  the  local 
gas  phase  sound  speed.  In  addition,  numerical  studies  indicate  that  for  D  >  D^j  the  gas 
velocity  relative  to  the  wave  head  is  locally  subsonic  at  the  non-ideal  strong  point,  while  the 
gas  velocity  relative  to  the  wave  head  is  locally  supersonic  at  the  non-ideal  weak  point. 
This  result  agrees  with  the  results  of  the  simple  one-phase  theory. 
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Shock  Discontinuity  Conditions 

A  shock  discontinuity  is  an  integral  part  of  a  two-phase  detonation.  As  in  one-phase 
ZND  theory  the  shock  wave  is  a  discontinuity  that  raises  the  pressure,  temperature,  and 
density  of  the  material,  initiating  significant  chemical  reaction.  In  the  context  of  the  one¬ 
dimensional  steady  model,  the  shock  wave  is  supported  by  the  chemical  energy  which  is 
released  by  the  reaction;  thus  the  process  is  self-sustaining. 

The  shock  conditions  are  determined  from  an  algebraic  analysis  and  provide  the  initial 
conditions  for  integrating  the  steady  equations  (5.1-4).  These  conditions  are  defined  by 
Equations  (5.1-15)  by  assuming  that  within  the  shock  wave,  reaction,  drag,  heat  transfer, 
and  compaction  have  no  effect.  Thus  through  the  shock  discontinuity,  differential 
equations  (5.1-4)  may  be  integrated  to  form  algebraic  relationships.  These  algebraic 
equations  admit  four  physical  solutions:  1)  the  ambient  state,  2)  shocked  gas,  unshocked 
solid,  3)  unshocked  gas,  shocked  solid,  and  4)  shocked  gas,  shocked  solid. 

This  model  ignores  the  effects  of  diffusive  momentum  and  energy  transport.  If 
included,  these  effects  would  define  a  shock  structure  of  finite  width.  Here  it  is  assumed 
that  the  length  scales  on  which  these  processes  are  important  are  much  smaller  than  the 
relaxation  scales  which  define  two-phase  detonation  structure.  To  the  author's  knowledge, 
this  assumption,  common  in  the  analysis  of  shocked  systems,  has  not  been  examined  either 
experimentally  or  theoretically  for  two-phase  reactive  systems. 

The  shock  conditions  are  given  below: 


-•s 

=  0 


(5.29) 

(5.30) 


(5.31) 
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(5.32) 
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Here  "s"  denotes  the  shocked  state  and  "0"  the  undisturbed  state.  Equations  (5.29-32) 
and  state  relations  (5.12, 13)  are  sufficient  to  calculate  the  shock  properties  for  phase  two. 
The  shock  state  for  the  solid  phase  is  independent  of  the  initial  porosity.  This  is  apparent 
from  Equation  (5.32),  which  says  that  the  porosity  does  not  change  through  the  shock 
discontinuity  (4)25  =  ^q),  and  from  Equations  (5.29-31)  where  it  is  seen  that  a  common 
factor  (|)2o  cancels  from  each  equation.  For  the  solid  phase  there  are  two  solutions  to 
Equations  (5.29-32):  the  inert  solution  and  the  shock  solution.  Exact  expressions  for  the 
solid  phase  shock  state  are: 
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In  Figure  5.4,  the  dimensional  solid  phase  shock  pressure  is  plotted  versus  the  shock  wave 
speed. 


D(m/s) 

Figure  5.4  Solid  Shock  Pressure  vs.  Shock  Wave  Speed 
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The  shock  properties  for  phase  one  are  inoiplied  by  the  mixture  equations  (5.5-7)  and 
state  relations  (5.9-10).  By  subtracting  the  solid  shock  equations  (5.29-31)  from  the 
respective  mixture  equations  (5.5-7),  one  obtains  gas  shock  jump  equations  which  are 
dependent  only  on  gas  phase  properties.  As  for  the  solid  phase,  in  these  equations  a 
conunon  factor  of  (|)io  cancels  from  each  equation.  Three  solutions  to  the  gas  shock  jump 
relations  exist:  the  inert  solution,  a  nonphysical  solution,  and  a  shock  solution.  In  the  limit 
as  7Ci3  (or  b)  approaches  zero,  the  nonphysical  prediction  of  gas  density  approaches  -l/7ti3 
and  is  therefore  rejected.  The  full  solution  to  the  non-ideal  shocked  gas  equations  are 
lengthy,  so  here  the  shocked  gas  solution  in  the  limit  of  an  ideal  gas  will  be  presented.  The 
full  non-ideal  shocked  solution  is  determined  by  solving  a  cubic  equation  described  in 
Appendix  E.  The  shocked  ideal  gas  state  (b  =  =  0)  is  described  by  the  following 

equations,  which  can  be  easily  rewritten  as  classical  ideal  gas  shock  relations  by  using  the 
definitions  of  the  dimensionless  parameters  to  write  these  equations  in  dimensional  form. 
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Figure  5.5  shows  a  plot  of  the  dimensional  gas  phase  shock  pressure  versus  the  shock 
wave  speed  for  the  non-ideal  gas. 
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D(m/s) 

Figure  5.5  Non-Ideal  Gas  Phase  Shock  Pressure  vs.  Shock  Wave  Speed 
Two-Phase  Detonation  Structure 

Before  studying  solutions  of  the  full  equations  (5.1-15),  a  simplified  model,  reduced 
to  two  differential  equations  is  considered.  These  equations  have  a  clear  geometrical 
interpretation  in  the  two-dimensional  phase  plane.  Results  from  this  model  will  be 
compared  to  those  of  the  full  model.  In  this  section  the  steps  necessary  to  reduce  Equations 
(5.1-15)  to  two  equations  will  be  described.  Next,  a  comparison  of  acceptable  detonation 
structure  predicted  by  the  two-equation  and  full  model  equations  is  given.  Finally,  an 
example  is  given  of  an  non-physical  solution  again  comparing  the  results  of  the  two- 
equation  model  with  those  of  the  full  model,  and  an  explanation  is  given  for  why  this 
solution  is  non-physical. 

The  steady  equations  (5.1-15)  are  simplified  significantly  when  heat  transfer  and 
compaction  effects  are  ignored.  This  corresponds  to  the  limit  Ttg  0  and  7C9  — ^  0.  From 
the  definition  of  the  dimensionless  parameter  713,  it  can  be  concluded  that  in  this  limit  the 
heat  transfer  in  the  reaction  zone,  roughly  h  /  D,  is  small  compared  to  the  thermal 
capacity  P2oCvi-  By  setting  7C9  to  zero,  it  is  assumed  that  compaction  effects  are  negligible; 
in  this  limit  Equation  (5.4)  holds  that  all  volume  change  is  due  to  chemical  reaction.  This  is 
achieved  mathematically  by  allowing  the  compaction  viscosity  |ip  to  approach  infinity. 

In  these  limits  it  is  possible  to  integrate  Equations  (5.1)  and  (5.3)  and  write  two 
autonomous  ordinary  differential  equations  in  two  unknowns,  solid  density  and  solid 
volume  fraction,  that  determine  the  system  completely.  All  other  thermodynamic  variables 
can  be  expressed  as  algebraic  functions  of  solid  density  and  volume  fraction.  With  the  two 
ordinary  differential  equations  it  is  easy  to  study  the  geometry  of  the  two-dimensional 
phase  space  in  the  p2-<|>2  phase  plane.  The  geometry  of  the  phase  plane  determines  whether 
a  detonation  structure  can  exist  in  theory. 
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To  derive  the  two-equation  model  requires  a  lengthy  algebraic  analysis.  Details  can  be 
found  in  Appendix  E.  To  summarize  the  process,  state  relations  (5.9,  10, 12, 13)  are  used 
to  eliminate  energy  and  temperamre  of  both  phases  in  all  remaining  equations.  Number 
conservation  (5.8)  is  used  to  eliminate  particle  radius  r  from  all  equations.  Mixture 
equations  (5.5-7)  are  used  to  write  gas  phase  properties  as  algebraic  functions  of  solid 
phase  properties.  In  uncoupling  the  mixture  equations,  a  complicated  cubic  equation  must 
be  solved.  One  root  corresponds  to  a  shocked  gas,  associated  with  what  is  known  in  one- 
phase  ZND  theory  as  the  strong  solution.  Another  root  corresponds  to  an  unshocked  gas, 
associated  with  the  weak  solution  in  one-phase  ZND  theory.  The  third  root  is  a  non¬ 
physical  consequence  of  the  virial  equation  of  state;  negative  gas  density,  temperature,  and 
pressure  are  predicted  with  this  root.  Substitution  of  these  results  into  Equations  (5.1-4) 
yields  four  ordinary  differential  equations  in  four  unknowns,  p2,  <1)2.  V2,  and  P2. 

When  the  limit  7C3  — >  0  and  JC9  0  is  considered,  combinations  of  two  of  these 
equations  can  be  integrated.  By  eliminating  the  gradient  of  volume  fraction  by  substituting 
Equation  (5.4)  into  (5.1),  a  homogeneous  equation  is  found  for  the  product  of  solid  density 
and  velocity.  When  integrated  this  gives  an  algebraic  relation  between  particle  density  and 
velocity.  The  solid  energy  equation  (5.3)  can  then  be  written  in  terms  of  a  homogeneous 
relation  involving  only  solid  pressure  and  density  by  using  the  integrated  mass  equation  to 
eliminate  velocity.  Initial  conditions  are  applied  corresponding  to  either  a  shocked  or 
unshocked  solid  state.  These  integrated  equations  allow  both  solid  velocity  and  pressure  to 
be  written  as  functions  of  solid  density.  The  integrated  equations  are  given  below. 


with 


(5.43) 

(5.44) 


shocked  solid 
unshocked  solid 


Here  K  is  a  constant  which  depends  on  whether  the  initial  solid  state  is  shocked  or 
unshocked.  With  these  results,  the  momentum  equation  (5.2)  can  be  used  to  determine  an 
explicit  equation  for  the  derivative  of  solid  density.  This  equation  along  with  the 
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compaction  equation  (5.4)  form  the  two-equation  model.  The  equations  which  govern  the 

structure  are  written  below  as 

<jpj  _  I-,) 

<14  g(Pj> 

(5.45) 

<*>2  _  IXp^. 
g(P2»  <l>2) 

(5.46) 

with  f,  g,  and  h  defined  as  follows 

f(p2.  <1>2)  =  V^1P2^? 

(5.47) 

g(p2.<t>2)  =  1) 

(5.48) 

h(p2.  <1>2)  =  ^iP2^2^r 

(5.49) 

These  equations  are  expressed  in  terms  of  the  functions  f,  g,  and  h,  which  are 
functions  of  p2  ‘1>2  is  seen  from  Equations  (5.45,  47)  that  the  solid  density 

changes  in  response  to  drag  effects,  embodied  in  the  terms  multiplying  the  drag  parameter 

7^2,  and  chemical  reaction  effects,  embodied  in  terms  multiplying  the  reaction  parameter  tcj. 

0 

Drag  terms  are  inherently  present  in  the  momentum  equation,  (5.2),  from  which  Equation 
(5.45)  is  derived.  Reaction  effects  arise  since  the  momentum  equation  (5.2)  predicts 
changes  in  momentum  due  to  changes  in  volume  fraction.  By  substituting  the  volume 
fraction  equation  (5.4)  into  the  momentum  equation,  reaction  effects  are  introduced. 
Effectively  then  the  momentum  equation  predicts  that  solid  density  changes  in  response  to 
drag  and  chemical  reaction  in  the  two-equation  model.  From  Equations  (5.46,  49)  it  is 
seen  that  volume  fraction  changes  are  predicted  only  in  response  to  chemical  reaction. 
Potential  equilibrium  states  exist  when  f  and  h  are  simultaneously  zero.  From  the 
functional  form  of  f  and  h,  it  is  seen  that  this  corresponds  to  a  state  where  density  changes 
due  to  drag  are  balanced  by  density  changes  due  to  reaction  and  where  volume  fraction 
changes  are  zero  due  to  complete  reaction. 
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When  g(P2,<l>2)  =  0,  and  f,  h  0,  infinite  gradients  are  predicted.  The  condition  g  =  0 
is  either  a  complete  reaction  or  solid  phase  sonic  condition  as  described  below.  Appendix 
E  shows  in  detail  how  for  the  two-equation  model,  the  solid  phase  is  sonic  when  Equation 
(5.51)  holds. 


r  =  0  (5.50) 

-1 

(5.51) 

When  either  Equations  (5.50)  or  (5.51)  hold,  forcing  g  to  zero,  it  is  seen  from  Equation 
(5.49)  that  h  is  simultaneously  zero. 

The  condition  g(P2>^)  =  0  leads  to  difficulties  regarding  the  division  by  zero.  The 
difficulties  in  the  continuation  of  the  solution  through  the  g  =  0  state  are  removed  by 
introducing  a  new  path  variable  z  and  considering  ^  as  an  independent  variable  defined  as 
follows 


=  m-O  (5.52) 

In  terms  of  the  new  independent  variable  z  Equations  (5.45,  46)  are  transformed  to  the 
following  equations 


dp. 

(5.53) 

d* 

(5.54) 

Equations  (5.53,  54)  are  autonomous  in  the  P2-<|>2  phase  plane.  Equation  (5.52)  may  be 
thought  of  as  an  auxiliary  relationship  to  determine  ^  once  the  structure  defined  by  the 
above  equations  is  determined.  Whether  Equations  (5.53,  54)  should  be  integrated 
forward  or  backward  in  z  is  a  relevant  question.  The  equations  should  be  integrated  so  that 
4  goes  from  0  to  -««.  From  Equation  (5.52)  it  is  seen  that  the  direction  of  change  of  ^  with 
respect  to  z  depends  on  whether  the  solid  phase  is  subsonic  or  supersonic.  If  the  initial 
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1 

j 

state  of  the  solid  is  unshocked,  the  solid  is  locally  supersonic,  g  >  0,  and  a  negative  dz  j 

corresponds  to  a  negative  d^.  If  the  initial  state  of  the  solid  is  shocked,  the  flow  is  locally 
subsonic,  g  <  0,  and  a  positive  dz  must  be  chosen  to  recover  a  negative  d^. 

In  the  context  of  this  reduced  model  there  are  several  requirements  for  an  admissible 
detonation  structure.  An  admissible  steady  structure  is  defined  by  an  integral  curve  which 

I 

begins  at  the  initial  point  in  the  P2-<1>2  plane  and  travels  in  that  plane  to  an  equilibrium 
position  where  f  and  h  are  simultaneously  zero.  This  point  is  defined  by  the  intersection  of 
the  curves  f  =  0  and  h  =  0.  In  addition  further  restrictions  are  placed  on  the  solution.  It  is 
required  that  the  gas  and  solid  thermodynamic  variables  density,  pressure,  and  temperature, 
are  always  positive  and  real.  Also  it  is  required  that  all  physical  variables  are  single-valued 
functions  of  the  position  variable  Based  on  these  restrictions  parametric  conditions  can 
be  obtained  for  admissibility  of  a  detonation  solution. 

The  conditions  under  which  thermodynamic  variables  become  either  negative  or 
imaginary  are  checked  numerically.  By  examining  a  few  limited  cases,  it  has  been  found 
that  there  are  regions  in  the  P2-<1>2  plane  where  gas  phase  pressure,  density,  and  temperature 
are  negative.  These  regions  are  bounded  by  curves  in  the  p2-<t>2  plane  where  gas  density, 
pressure,  and  temperature  are  zero.  In  solving  the  cubic  equation  for  the  gas  phase 
propenies,  imaginary  gas  phase  quantities  are  sometimes  predicted.  It  has  been  found 
numerically  that  the  border  of  the  imaginary  region  corresponds  to  a  sonic  condition  in  the 
gas  phase. 

The  geometry  of  the  f  =  0  and  h  =  0  curves  is  critical  in  determining  the  integral  curve 
which  defines  the  steady  state  solution.  Depending  on  the  relative  orientation  of  these 
curves  and  the  initial  state,  many  classes  of  solutions,  each  with  a  distinct  character,  are 
available.  Some  solutions  reach  an  equilibrium  state,  defined  at  the  intersection  of  the  f  =  0 
and  h  =  0  curves.  The  structure  of  the  steady  detonation  solution  is  strongly  influenced  by 
the  nature  of  the  equilibrium  point,  which  can  be  classified  as  a  source,  sink,  saddle,  or 
spiral.  For  example,  if  the  equilibrium  point  to  which  the  integral  curve  is  drawn  is  a  sink, 
then  a  continuum  of  wave  speeds  are  found  for  which  steady  detonations  are  allowed.  If 
t.he  equilibrium  point  is  a  saddle,  there  is  only  one  wave  speed  which  will  bring  the  integral 
curve  to  the  equilibrium  position.  For  some  wave  speeds  the  orientation  of  the  f  =  0  and  h 
=  0  curves  prevents  solutions  from  reaching  an  equilibrium  state;  these  solutions  cannot  be 
classified  as  steady  solutions.  Among  these  types  of  solutions  are  those  that  pass  through  a 
solid  sonic  point  and  become  physically  unacceptable,  multivalued  functions  of  distance. 

Figure  5.6  shows  sketches  of  phase  planes  for  two  classes  of  solutions,  one  an  acceptable 
detonation  structure,  the  other  a  nonphysical  solution. 
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<l>2  Acceptable  Detonation 
A  Phase  Plane  Staicture 


Nonphysical  Solution 


Equilibrium 

Point 

g=h=f=0 


Sonic  Line 
g  =  h  =  0 


\  '2 
Complete  Reaction 

Line  g  =  h  =  0 


f  =  0 


Sonic  Line 
g  =  h  =  0 


Equilibrium  ^  . 

PqJjjj  Complete  ReacUon 

g  =  h  =  f=0  Lineg  =  h  =  0 


Figure  5.6  Phase  Plane  Sketches  of  Physical  and  Nonphysical  Solutions 


Each  sketch  shows  the  separatrix  lines  f  =  0  and  g  =  h  =  0.  The  equilibrium  position 
is  at  the  intersection  of  these  curves.  Each  curve  shows  a  solid  phase  sonic  line,  g  =  h  =  0, 
forbidden  regions  in  which  gas  phase  properties  are  not  physical,  and  integral  curves  which 
originate  from  the  initial  condition.  For  the  acceptable  structure  the  integral  curve  travels 
from  the  initial  state  to  the  equilibrium  position.  By  changing  the  flow  conditions,  the 
topology  of  this  phase  plane  is  altered,  shown  in  the  adjacent  sketch.  In  this  sketch,  the 
integral  curve  is  driven  through  the  solid  sonic  line  and  is  incapable  of  reaching  the 
equilibrium  point.  As  explained  below,  past  the  solid  sonic  line,  the  solution  is  double¬ 
valued  and  therefore  not  physical. 

Thermodynamic  variables  become  double-valued  functions  of  distance  when  a  solid 
sonic  condition  (g  =  0)  is  reached  at  a  non-equilibrium  point  in  the  phase  plane  (f  ^  0). 
From  Equation  (5.52),  it  is  seen  that  the  direction  of  change  of  %  with  respect  to  z  changes 
when  the  solution  passes  through  a  solid  sonic  point.  Thus  which  starts  at  zero  and 
moves  towards  -<»  as  reaction  progresses,  changes  direction  and  moves  towards  +<»  at  a 
critical  point  when  a  solid  sonic  condition  is  reached.  Through  this  point  Equations 
(5.53,  54)  predict  a  continuous  variation  of  density  and  volume  fraction.  At  the  solid  sonic 
point  the  derivatives  of  pj  and  <J>2  with  respect  to  z  are  finite,  and  the  derivatives  with 
respect  to  ^  are  infinite.  At  any  given  location  ^  two  values  of  each 

thermodynamic  variable  will  be  predicted.  This  is  physically  unacceptable. 
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Results  analogous  to  one-phase  ZND  theory  can  be  obtained  with  the  two-equation 
model.  For  the  input  conditions  of  Table  I,  with  the  heat  transfer  coefficient  h  =  0  and 
compaction  viscosity  Pc  — > «»,  an  initial  porosity  greater  than  0.19,  and  an  initially  shocked 
gas  and  unshocked  solid,  a  CJ  structure  can  be  defmed.  In  these  limits  there  is  no  heat 
transfer  or  volume  change  due  to  pore  collapse.  The  CJ  wave  speed  is  determined  from 
solving  the  earlier-described  equation  set  (5.21-23).  Wave  speeds  less  than  the  CJ  speed 
are  rejected  because  imaginary  gas  phase  quantities  are  predicted  near  the  complete  reaction 
end  state.  Wave  speeds  greater  than  CJ  are  admitted  by  this  model  and  correspond  to  the 
strong  ZND  solution.  Such  a  wave  leaves  the  gas  at  a  velocity  which  is  subsonic  relative  to 
the  wave  front.  As  in  ZND  theory,  piston  support  is  required  to  prevent  rarefaction  waves 
from  damping  the  reaction  zone  structure.  For  the  CJ  wave,  the  final  velocity  is  sonic  and 
no  piston  support  is  necessary  to  support  the  wave. 

The  solution  is  driven  to  a  sink  in  the  p2-<}>2  plane.  To  show  this  point  is  a  sink,  one 
first  finds  the  equilibrium  point  by  solving  the  algebraic  problem  f(p2,  4)2)  =  0,  h(p2,  (1)2)  = 
0.  The  differential  equations  (4.53,  54)  are  then  linearized  about  this  equilibrium  point. 
These  linear  differential  equations  can  be  solved  exactly  to  determine  the  behavior  of  any 
integral  curve  which  approaches  the  equilibrium  point.  In  this  study,  for  an  shtxked  sohd 
and  shocked  or  unshocked  gas,  it  was  found  that  all  integral  curves  in  the  neighborhood  of 
the  equilibrium  point  were  attracted  to  the  equilbrium  point;  in  the  terminology  of  ordinary 
differential  equation  theory,  that  point  is  classified  as  a  sink. 

The  ordinary  differential  equations  of  the  two-equation  model  and  full,  four-equation 
model  were  solved  numerically.  Integration  was  performed  using  the  IMSL  subroutine 
DVERK,  a  fifth  and  sixth  order  Runge-Kutta  routine,  on  the  UIUC  Cyber  175.  Step  sizes 
were  chosen  such  that  none  of  the  fundamental  variables,  P2,  <1>2>  ^2.  and  P2,  changed  by 
more  than  5%  in  value  in  any  given  integration  step.  Typically  about  two  hundred 
integration  steps  v/ere  sufficient  to  describe  the  reaction  zone.  A  typical  integration  took 
twenty  seconds  to  complete. 

For  an  initial  solid  volume  fraction  of  0.70,  Figure  5.7  shows  a  plot  of  the  phase  plane 
for  a  CJ  wave  speed  of  7369  m/s.  This  curve  shows  the  sonic  line  (g  =  h  =  0)  on  P2  = 
1.35,  the  complete  reaction  line  (g  =  h  =  0)  on  <()2  =  0  and  the  f  =  0  line.  It  is  seen  from  this 
curve  that  the  only  equilibrium  point  is  at  -  (1-04,  0).  The  vector  field 

superimposed  on  this  figure,  defined  by  Equations  (5.53,  54),  shows  this  point  is  a  sink 
which  is  confirmed  by  a  local  linear  analysis  near  the  equilibrium  point.  The  integral  curve 
connecting  the  initial  state  to  the  equilibrium  point  is  also  plotted  on  this  figure.  This  curve 
is  obtained  by  numerical  integration  of  Equations  (5.53,  54).  This  integral  curve  moves  in 
a  direction  defined  by  the  vector  field  of  the  phase  plane.  Curves  of  zero  gas  phase 
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pressure  are  plotted  in  this  figure  along  with  the  curve  defining  the  boundary  between  pure 
real  and  imaginary  gas  phase  quantities.  The  gas  velocity  is  locally  sonic  (Mj^  =  1)  on  the 
boundary  of  the  region  where  imaginary  gas  phase  properties  exist.  This  indicates  that  the 
solution  is  non-physical  if  the  gas  passes  through  a  sonic  condition  at  a  point  of  incomplete 
reaction. 

When  the  full  model  equations  are  considered,  general  results  from  the  two-equation 
model  are  retained.  It  is  more  difficult  to  interpret  these  results  as  the  phase  space  is  four¬ 
dimensional.  With  a  given  set  of  initial  conditions,  the  gas  phase  CJ  end  state  is  the  same 
whether  the  two-equation  or  four-equation  model  is  used.  The  solid  phase  end  state  and 
details  of  the  reaction  zone  structure  do  depend  on  which  model  equations  are  used.  Plots 
of  predicted  detonation  structure  are  shown  in  Figure  5.8,  which  plots  solid  and  gas 
density,  lab  velocity  u,  pressure,  temperature,  Mach  number,  particle  radius,  and  solid 
volume  fraction  versus  distance  Also  plotted  on  this  figure  are  results  from  the  two- 
equation  model.  It  is  seen  that  both  models  predict  results  of  the  same  order  of  magnitude. 
Gas  phase  quantities  are  nearly  identical  for  both  models.  While  there  are  small  differences 
in  solid  phase  predictions,  these  results  are  remarkable  as  there  is  no  real  basis  to  assume 
the  the  limits  taken  are  appropriate  for  this  class  of  models.  These  results  show  that 
material  compaction  and  heat  transfer  are  not  important  mechanisms  in  determining  two- 
phase  detonation  structure  and  that  there  is  justification  in  using  the  two-equation  model  as 
a  tool  for  understanding  the  full  model  equations.  A  comparison  of  some  results  of  the  two 
models  is  given  in  Table  H. 

In  Figure  5.8  it  is  seen  that  the  gas  phase  is  shocked  while  the  solid  phase  is 
unshocked.  It  is  noted  from  Figure  5.8c  that  the  gas  pressure  continues  to  rise  past  the 
initial  shock  gas  pressure,  in  contrast  to  the  results  of  the  simple  one-phase  ZND  theory, 
which  predicts  the  pressure  to  be  a  maximum  at  the  shock  state.  From  this  maximum, 
known  as  the  "Von  Neumann  spike,"  the  pressure  decreases  to  the  equilibrium  CJ 
pressure.  It  should  also  be  noted  that  the  high  gas  phase  temperature  (~  10,000  K) 
indicates  that  ionization,  dissociation,  and  radiative  heat  transfer  could  be  important 
mechanisms  in  the  reaction  zone.  These  effects  have  not  been  considered  but  could  be 
incorporated  into  future  work. 

Non-physical  solutions  are  now  considered.  Such  solutions  exist  below  a  critical 
value  of  initial  solid  volume  fraction.  The  critical  point  is  shown  in  Figure  5.9,  which  plots 
CJ  wave  speed  versus  initial  bulk  density  (p^  =  Pio<)>io  Pzo^bo)-  figure  also 
compares  predictions  of  this  model  with  those  of  the  unsteady  model  of  Butler  and  Krier 
[  1  ]  and  those  of  the  equilibrium  thermochemistry  code  TIGER  given  in  Ref.  1 .  The  feature 
of  a  critical  initial  bulk  density  has  not  been  identified  by  other  models. 


Figure  5.8c  Gas  and  Solid  Pressure  Two-Phase  CJ  Detonation  Structure 


Figure  5.8d  Gas  and  Solid  Temperature  Two-Phase  CJ  Detonation  Structure 
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Table  H 

COMPARISON  OF  TWO  AND  FOUR  EQUATION  MODEL  PREDICTIONS  FOR 

CJ  Waves  with  and  without  Leading  Gas  phase  Shock 

Leading  Gas  Shock  No  Leading  Gas  Shock 


Two-Equation 

Full 

Two-Equation 

EnU 

Initial  Bulk  Density 

1,333  kg/m3 

1,333  kg/m3 

1,333  kg/m3 

1,333  kg/m3 

Reaction  Zone  Length 

13.(X)  mm 

12.89  mm 

62.1  mm 

61.7  mm 

CJ  Wave  Speed 

7,369  m/s 

7,369  m/s 

7,369  m/s 

7,369  m/s 

CJ  Pressure 

19.4  GPa 

19.4  GPa 

19.4  GPa 

19.4  GPa 

CJ  Density 

1,821  kg/m3 

1,821  kg/m3 

1,821  kg/m3 

1,821  kg/m3 

CJ  Temperature 

4,176  K 

4.176  K 

4,174  K 

4,174  K 

CJ  Gas  Velocity 

1,976  m/s 

1.976  m/s 

1,974  m/s 

1,974  m/s 

(CJ  Gas  Mach  Number)^ 

1 

1 

1 

1 

Maximum  Gas  Temperature 

11,119  K 

11.108  K 

4,174  K 

4,174  K 

Final  Solid  Pressure 

0.716  GPa 

0.636  GPa 

0.716  GPa 

0.636  GPa 

Final  Solid  Density 

1,973  kg/m3 

1.962  kg/m3 

1,973  kg/m3 

1,962  kg/m3 

Final  Solid  Temperature 

349  K 

344K 

349  K 

344  K 

Final  Solid  Velocity 

272  m/s 

429  m/s 

272  m/s 

428  m/s 

(Final  Solid  Mach  Number)^  4.82 

4.67 

4.82 

4.67 
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Figure  5.9  CJ  Wave  Speed  vs.  Initial  Bulk  Density 


For  a  value  of  initial  solid  volume  fraction  of  0.20,  very  near  the  critical  bulk  density, 
an  acceptable  detonation  structure  is  obtained.  A  phase  portrait,  vector  map,  and  integral 
curve  is  shown  in  Figure  5.10.  The  figure  resembles  Figure  5.6,  but  the  curves  have  all 
been  skewed.  Note  that  the  integral  curve  nearly  reaches  the  sonic  state  before  turning 
around  and  travelling  to  the  complete  reaction  end  state. 

For  an  initial  solid  volume  fraction  of  0.15,  a  non-physical  solution  is  obtained  for  a 
CJ  wave  speed.  The  two-equation  model's  phase  plane  is  shown  in  Figure  5.11.  The 
integral  curve  in  this  plane  passes  through  the  solid  sonic  line  at  a  non-equilibrium  point 
causing  the  solution  to  become  double-valued.  A  plot  of  the  solid  phase  Mach  number  is 
shown  in  Figure  5.12  for  both  the  two  and  four  equation  models.  Again  both  models 
predict  nearly  identical  results.  It  is  seen  from  Figure  5.12  that  infinite  gradients  with 
respect  to  ^  are  predicted  precisely  at  the  point  where  the  solid  phase  reaches  a  sonic 
velocity  (M2^  =  1). 
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Figure  5.12  Solid  Phase  Mach  Number  for  Nonphysical  Solution 

Solutions  with  no  leading  shock  in  either  the  gas  or  solid  phase  are  also  admitted  by 
this  model.  Figure  5.13  shows  the  phase  portrait,  vector  map  and  integral  curve  for  a  CJ 
wave  with  no  leading  gas  or  solid  shock  propagating  through  a  mixture  with  an  initial  solid 
volume  fraction  of  0.70.  Again,  the  equilibrium  point  is  a  sink.  As  summarized  in  Table 
II,  the  main  difference  between  this  case  and  the  case  with  the  leading  gas  phase  shock  is 
that  the  reaction  zone  is  much  longer  (62  mm  vs.  13  mm)  for  no  leading  shock  in  the  gas 
phase.  Again  both  two  and  four  equation  models  predict  similar  results.  The  CJ  gas  phase 
end  state  is  identical  regardless  of  whether  the  initial  gas  state  is  shocked  or  unshocked,  or 
whether  the  two  or  four  equation  model  is  used.  This  is  because  the  complete  reaction  CJ 
state  is  independent  of  the  structure  of  the  detonation.  Small  differences  in  the  CJ 
temperatures  and  gas  velocities  can  be  attributed  to  numerical  roundoff  errors  as  the  CJ 
state  is  extremely  sensitive  to  the  CJ  wave  speed.  In  general  the  solid  end  state  can  vary  for 
each  state  presented  in  Table  H.  It  is  noted  that  the  solid  phase  end  state  predicted  by  the 
two-equation  model  is  nearly  the  same  for  both  the  un  shocked  and  shocked  gas  as  is  the 
solid  phase  end  state  for  the  full  equation  model. 
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For  wave  speeds  greater  than  CJ,  strong  and  weak  waves  can  be  predicted.  For  an 
initial  solid  volume  fraction  of  0.70  and  a  wave  speed  of  8,000  ir  /s  (which  is  greater  than 
the  CJ  wave  speed  of  7369  m/s)  Figures  5.14  and  5.15  show  plots  of  the  two-equation 
model's  phase  portraits  for  the  strong  (initially  shocked  gas)  and  weak  (initially  unshocked 
gas)  case.  In  each  case  the  solid  is  initially  unshocked.  The  equilibrium  points  are  sinks  in 
both  cases.  The  results  of  these  calculations  for  both  two  and  four  equation  models  are 
summarized  in  Table  III.  For  the  strong  case  the  reaction  zone  is  shorter  than  for  the 
corresponding  CJ  wave  with  a  leading  gas  phase  shock.  For  the  weak  case  the  reaction 
zone  is  longer  than  for  the  corresponding  CJ  wave  without  a  leading  gas  phase  shock. 
Again  two  and  four  equation  models  predict  similar  results. 

This  study  predicts  a  continuum  of  two-phase  detonation  wave  speeds  as  a  function  of 
piston  velocity.  CJ  wave  speed  is  plotted  as  a  function  of  piston  velocity  in  Figure  5.16. 
For  wave  speeds  greater  than  CJ,  piston  support  is  required  to  support  the  wave.  The  CJ 
wave  can  propagate  with  or  without  piston  support  as  the  complete  reaction  point  is  a  gas 
phase  sonic  point. 


Figure  5.16  Two-Phase  Detonation  Wave  Speed  vs.  Piston  Velocity 

For  piston  velocities  below  CJ,  a  continuum  of  weak  waves  are  predicted.  The 
implications  of  this  are  unclear.  As  the  complete  reaction  point  is  supersonic,  the  piston 
support  is  not  necessary.  This  suggests  that  the  solution  may  not  be  unique.  Simple  one- 
phase  ZND  theory  also  predicts  a  continuum  of  weak  waves.  Fickett  and  Davis  [40] 
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Table  HI 


COMPARISON  OF  TWO  AND  FOUR  EQUATION  MODEL  PREDICTIONS  FOR 

Strong  and  weak  detonations,  d  =  8,000  m/s 

Leading  Gas  Shock  CSirong)  No  Leading  Gas  Shock  (Weak) 


Two-Eauation 

Full 

Two-Equation 

Full 

Initial  Bulk  Density 

1,333  kg/m^ 

1,333  kg/m3 

1,333  kg/m3 

1,333  kg/m3 

Reaction  Zone  Length 

10.2  mm 

10.1  mm 

71.4  mm 

70.9  mm 

Wave  Speed 

8,000  m/s 

8,000  m/s 

8,000  m/s 

8,000  m/s 

Final  Gas  Pressure 

32.3  GPa 

32.3  GPa 

13.8  GPa 

13.8  GPa 

Final  Gas  Density 

2,145  kg/m3 

2,145  kg/m3 

1,590  kg/m3 

1 ,590  kg/m3 

Final  Gas  Temperature 

5,274  K 

5,274  K 

3,710  K 

3,710  K 

Final  Gas  Velocity 

3,029  m/s 

3,029  m/s 

1,291  m/s 

1,291  m/s 

(CJ  Gas  Mach  Number)^ 

0.567 

0.567 

1.99 

1.99 

Maximum  Gas  Temperature 

12,526  K 

12,514  K 

3,710  K 

3,710  K 

Final  Solid  Pressure 

0.744  GPa 

0.676  GPa 

.657  GPa 

.542  GPa 

Final  Solid  Density 

1,975  kg/m3 

1,966  kg/m3 

1,967  kg/m3 

1,953  kg/m3 

Final  Solid  Temperature 

351  K 

346  K 

345  K 

337  K 

Final  Solid  Velocity 

306  m/s 

484  m/s 

273  m/s 

432  m/s 

(Final  Solid  Mach  Number)'' 

-  5.63 

5.44 

5.78 

5.66 

I 
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discuss  this  issue  for  one-phase  theory.  Though  this  issue  is  still  not  settled  for  the  one- 
phase  model,  some  have  suggested  that  the  weak  waves  may  be  ruled  out  as  unphysical 
because  of  a  lack  of  an  initiation  mechanism.  Pickett  and  Davis  show  results  of  more 
complicated  one-phase  models  which  indicate  that  a  unique  weak  wave  speed  exists  when 
such  mechanisms  as  diffusive  heat  and  momentum  transfer  are  taken  into  account.  A 
similar  result  may  hold  for  two-phase  detonations. 
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VI.  Conclusions  and  recommendations 


Compaction  Waves 

The  piston-impact  problem  for  a  compressible  porous  solid  has  been  solved  in  the 
context  of  a  steady  two-phase  model  neglecting  gas  phase  effects.  With  this  model,  it  is 
possible  to  obtain  an  exact  solution  for  the  compaction  wave  speed,  final  porosity,  and  final 
pressure.  The  degree  of  accuracy  of  the  predictions  can  be  attributed  to  the  ad  hoc 
estimates  for  the  non-ideal  solid  parameter  and  the  assumed  form  of  the  static  pore  collapse 
function,  f.  Within  the  framework  of  this  model  it  is  possible  to  understand  the  general 
features  of  a  compaction  wave.  Two  classes  of  compaction  waves  have  been  identified, 
subsonic  waves  with  no  leading  shock,  and  supersonic  waves  with  a  leading  shock.  It  is 
predicted  that  the  magnitude  of  the  supporting  piston  velocity  determines  which  class  of 
wave  exists,  with  low  piston  velocities  resulting  in  a  subsonic  structure  and  high  piston 
velocities  resulting  in  a  supersonic  structure. 

A  compaction  wave  with  structure  has  been  predicted  because  a  dynamic  pore  collapse 
equation  has  been  used.  As  summarized  by  Kooker  [62],  many  compaction  wave  models 
do  not  consider  dynamic  pore  collapse;  rather  they  enforce  static  pore  collapse  (P  =f) 
throughout  the  flow  field.  In  zero  gas  density  limit,  such  an  assumption  results  in  a 
compaction  wave  without  structure.  The  pressure  discontinuously  adjusts  to  a  static 
equilibrium  value.  However,  it  is  not  established  whether  two-phase  models  with  static 
pore  collapse  are  hyperbolic,  a  necessary  condition  if  discontinuities  are  to  be  admitted  and 
for  a  well-posed  initial  value  problem.  For  two-phase  models  assuming  pressure 
equilibrium  between  phases  but  not  incorporating  quasi-static  compaction  data, 
Lyezkowski,  et.  al.  [53]  have  identified  regimes  in  which  unsteady  two-phase  equations 
are  not  hyperbolic. 

There  are  many  ways  to  extend  the  compaction  wave  study.  By  including  the  effects 
of  the  gas  phase,  it  should  be  possible  to  determine  how  the  gas  phase's  presence  modifies 
the  compaction  wave  structure.  By  including  the  effect  of  particle  size  in  f,  it  should  be 
possible  to  model  the  experiments  of  Elban,  et  al.  [63]  which  show  that  the  static  pore 
collapse  stress  level  is  a  function  of  both  volume  fraction  and  particle  size.  By  considering 
the  solid  to  be  composed  of  particles,  it  may  be  possible  to  model  the  effect  of  particle 
breakup  on  the  results  when  f  is  assumed  to  be  a  function  of  particle  size. 


Detonation  Wa\)es 


It  is  thought  that  the  most  important  contribution  of  this  study  is  that  existence 
conditions  have  been  predicted  for  a  steady,  one-dimensional,  two-phase  detonation  in  a 
granular  material.  The  available  detonation  solutions  are  restricted  by  both  algebraic 
equilibrium  end  state  analysis  and  by  an  analysis  of  the  structure  of  the  steady  wave. 

Though  gas  phase  end  state  analysis  has  been  performed  by  many  others,  it  is  believed 
that  the  work  here  clarifies  this  analysis  by  finding  simple  analogies  between  one-phase  CJ 
conditions  and  two-phase  CJ  conditions  along  with  simple  corrections  for  non-ideal  gas 
phase  effects.  These  simple  two-phase  conditions  are  analogous  to,  but  not  identical  to,  the 
one-phase  CJ  condition  and  cannot  be  obtained  a  priori  from  the  one-phase  model.  The 
similarity  in  results  is  due  to  the  similarities  which  exist  between  the  one-phase 
conservation  equations  and  two-phase  conservation  equations.  The  common  notion  that 
one-phase  CJ  results  can  be  directly  applied  to  two-phase  systems  is  disproved  by  this 
work. 

The  variation  of  CJ  properties  with  initial  bulk  density  reported  here  accurately 
matches  the  TIGER  predictions  for  a  single  set  of  gas  phase  state  parameters.  Thus  it  is  not 
necessary  to  vary  the  gas  phase  state  equation  parameters  as  initial  bulk  density  changes  to 
match  the  TIGER  predictions  as  done  by  other  researchers.  In  Ref.  47  a  virial  equation  of 
state  identical  in  form  to  the  gas  state  equation  of  this  study  was  used.  In  that  study  as  the 
initial  bulk  density  varied,  the  value  of  b  was  varied  within  the  range  from  0.00361  m^/kg 
to  0.00486  m^/kg  in  order  to  match  the  TIGER  predictions.  As  shown  in  Figure  5.3  of  the 
present  study  the  CJ  properties  are  very  sensitive  to  changes  in  b  on  the  order  of  those 
studied  in  Ref.  47.  In  Ref.  2  a  JWL  gas  state  equation  is  used,  and  it  is  reported  that  CJ 
data  is  adequately  reproduced  when  the  constants  are  allowed  to  vary  with  the  initial  bulk 
density.  It  is  believed  that  the  approach  of  the  present  study  in  determining  CJ  properties 
has  the  advantage  over  the  approach  taken  in  Refs.  47  and  2.  Though  all  the  studies  fix  gas 
state  equation  parameters  so  that  CJ  predictions  or  data  is  matched,  a  single  set  of 
parameters  is  used  only  in  the  present  study. 

An  analysis  of  the  structure  of  a  two-phase  detonation  wave  has  further  restricted  the 
class  of  available  steady  solutions.  The  structure  analysis  has  shown  that  below  a  critical 
initial  bulk  density  no  steady  solution  can  exist  when  the  solid  particles  reach  a  sonic  state. 
The  mathematical  consequence  of  this  is  that  the  solution  is  becomes  a  double-valued 
function  of  distance,  a  physically  unacceptable  result.  This  particular  result  and  the  general 
technique  of  using  structure  analysis  to  limit  the  available  solutions  is  new  to  two-phase 
detonation  theory. 
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As  a  result  of  this  study  it  is  possible  to  predict  the  features  of  a  steady  two-phase 
detonation  structure.  It  has  been  shown  that  when  a  leading  shock  wave  exists  in  the  gas 
phase  and  the  solid  is  unshocked,  that  two-phase  equivalents  to  the  one-phase  ZND  strong 
and  CJ  solutions  are  predicted.  As  in  ZND  theory,  the  two-phase  theory  predicts  that 
piston  support  is  required  for  the  strong  solution  to  exist,  and  that  a  two-phase  CJ 
detonation  can  propagate  wi^h  or  without  piston  support.  It  has  also  been  shown  that  when 
both  the  gas  and  solid  phases  are  unshocked,  that  the  model  equations  yield  two-phase 
equivalents  of  weak  and  CJ  solutions.  These  types  of  solutions  are  also  found  using  the 
simple  one-phase  ZND  theory  but  are  comriKsnly  dismissed  because  it  is  thought  there  is  no 
mechanism  to  initiate  reaction.  The  model  yields  such  solutions  because  the  functional 
form  of  the  combustion  model  allows  a  small  amount  of  reaction  to  occur  even  at  ambient 
conditions.  The  model  allows  the  small  heat  released  by  the  reaction  to  accumulate  and 
cause  a  thermal  explosion  after  an  induction  time. 

This  work  has  clarified  the  role  of  shock  jumps  in  two-phase  detonation  theory.  No 
previous  work  on  two-phase  detonation  theory  has  considered  the  four  possible  states 
admitted  by  the  shock  jump  conditions.  This  study  has  shown  that  two-phase  detonation 
structure  is  possible  when  the  gas  phase  is  shocked  or  unshocked  and  the  solid  is 
unshocked.  The  possibility  of  a  two-phase  detonation  with  a  shocked  solid  has  not  been 
ruled  out;  an  example  of  such  a  detonation  has  not  been  found  yet.  This  study  does  not 
consider  how  the  structure  of  an  unshocked  solid  and  shocked  gas  can  arise.  To  show 
how  this  could  occur  would  require  an  unsteady  analysis  which  is  beyond  the  scope  of  this 
study. 

To  speculate  on  how  such  a  scenario  could  develop,  on  could  imagine  a  slow, 
unconfined  burning  of  reactive  particles.  If  the  system  were  suddenly  confined,  a  local 
region  of  high  gas  pressure  could  develop  which  could  give  rise  to  a  propagating  shock 
wave  in  the  gas  but  not  the  solid.  It  should  also  be  said  that  the  idea  of  shocked  gas  and 
unshocked  solid  is  common  in  the  literature  of  shock  waves  in  dusty  gases.  A  standard 
assumption  is  that  there  is  a  shock  wave  in  the  gas  but  that  the  solid  particles  are 
incompressible,  thus  unshocked.  Rudinger  [64]  provides  an  example  of  suc..  a  model. 

This  study  has  for  the  first  time  unambiguously  identified  a  finite-valued  gas  and  solid 
complete  reaction  end  state.  Though  others  have  discussed  the  gas  phase  complete  reaction 
end  state,  the  solid  end  state  has  never  been  considered.  In  each  of  the  physical  detonation 
solutions  presented  here  the  final  values  of  both  the  solid  and  gas  can  be  precisely  stated. 
In  all  cases,  the  complete  reaction  end  state  analysis  allows  the  final  gas  phase  properties  to 
be  determined.  For  the  two-equation  model,  the  final  solid  properties  can  be  determined  by 
an  algebraic  analysis  without  regard  to  the  detonation  structure. 
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The  complete  reaction  singularity  which  exists  due  to  the  1/r  terms  in  the  governing 
equations  leads  one  to  question  whether  unbounded  properties  are  predicted  at  complete 
reaction.  Previous  studies  have  neglected  this  question.  Here  it  has  been  shown  that  a 
two-phase  detonation  can  be  predicted  when  p.-oper  account  is  taken  for  the  complete 
reaction  singularity. 

This  study  has  also  identified  for  the  first  time  the  importance  of  sonic  singularities  in 
two-phase  detonation  systems.  It  has  been  shown  that  in  general  if  a  sonic  condition  is 
reached  in  the  solid  phase,  that  double-valued  properties  are  predicted,  and  that  if  a  sonic 
condition  is  reached  in  the  gas  phase  at  a  point  of  incomplete  reaction,  that  imaginary  gas 
phase  properties  are  predicted.  The  sonic  conditions  are  particular  for  each  phase  and  have 
no  relation  to  the  mixture  sound  speed. 

Techniques  which  are  new  in  the  two-phase  detonation  modeling  field  have  been  used 
to  simplify  the  governing  equations.  An  algebraic  method  for  uncoupling  the  mixture 
mass,  momentum,  and  energy  equations  to  solve  for  gas  phase  variables  in  terms  of  solid 
phase  variables  has  been  developed.  It  has  been  shown  that  the  equations  can  be  reduced 
to  a  set  of  four  uncoupled  ordinary  differential  equations  in  four  unknowns  and  how  in  the 
limit  of  zero  heat  transfer  and  compaction  these  equations  reduce  to  two  ordinary 
differential  equations.  The  two-equation  model  makes  it  possible  to  exploit  the  simple  two- 
dimensional  phase  plane  to  gain  understanding  of  the  complete  model.  Similarity  of  the 
results  of  the  two  and  four  equation  models  suggests  that  heat  transfer  and  compaction  are 
not  important  mechanisms  in  determining  two-phase  detonation  structure. 

Much  work  remains  to  be  done  in  two-phase  detonation  theory.  It  is  highly  likely  that 
other  classes  of  steady  detonations  can  be  predicted  which  have  not  been  studied  here.  The 
complexity  of  the  model  equations  makes  this  search  a  trial  and  error  process.  However, 
one  can  envision  several  different  detonation  scenarios  by  making  minor  adjustments  in  the 
relative  positions  of  the  separatrices  in  the  two-dimensional  phase  plane. 

Two-phase  steady  detonation  results  can  be  effectively  used  in  the  unsteady  two-phase 
DDT  problem.  Predictions  of  any  unsteady  model  would  be  strengthened  by  comparing 
them  to  the  predictions  of  a  steady  model.  Unsteady  model  results  can  be  used  to  verify 
that  the  unsupported  two-phase  detonation  wave  is  a  CJ  wave.  This  would  simply  require 
an  examination  of  the  two-phase  end  state  conditions. 

Reaction  zone  lengths  predicted  by  the  steady  model  must  match  those  predicted  by  the 
unsteady  model.  This  however  raises  an  important  question  regarding  numerical 
resolution.  This  study  predicts  reaction  zone  lengths  of  the  order  of  10  mm.  Unsteady 
two-phase  models  now  use  a  cell  size  on  the  order  of  1  mm.  It  is  highly  unlikely  that  with 
the  ratio  of  cell  size  to  reaction  zone  length  so  high  that  one  could  use  unsteady  results  to 
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distinguish  features  of  the  reaction  zone  identified  by  steady  analysis,  in  particular,  shock 
waves.  The  results  arc  smeared  by  artificial  viscosity  and  lack  of  an  adequate  number  of 
cells.  Thus  the  results  of  this  study  suggest  that  a  cell  size  on  the  order  of  0.01-0. 1  mm  be 
employed  in  unsteady  calculations.  Cell  sizes  of  this  magnitude  present  a  dilemma. 
Typical  particle  sizes  for  detonation  applications  range  from  0.1-1  mm.  One  assumption  of 
continuum  modeling  of  granular  materials  is  a  large  number  of  particles  exist  in  any 
averaging  volume.  If  cell  sizes  of  the  order  of  0.01-0.1  mm  are  employed,  as  the  results 
suggest  is  necessary,  then  the  continuum  assumptions  may  not  be  valid. 

The  results  of  two-phase  steady  theory  can  be  used  as  the  basis  for  further  studies.  At 
this  time,  the  stability  of  two-phase  detonations  has  yet  to  be  investigated.  Also 
multidimensional  two-phase  theory  is  undeveloped.  It  may  be  possible  to  obtain  a 
relationship  to  determine  the  critical  diameter  of  a  cylinder  containing  a  two-phase  explosive 
much  in  the  same  way  these  relations  have  been  developed  for  one-phase  materials  [65]. 
Finally,  it  should  be  possible  to  use  the  method  of  characteristics  to  study  the  unsteady 
two-phase  problem  in  a  new  way  which  has  the  potential  to  provide  more  understanding  of 
what  processes  actually  cause  a  two-phase  detonation. 
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Appendix  A.  Characteristic  Form  of  Governing  equations 


This  appendix  will  identify  the  characteristic  directions  and  characteristic  form  of 
Equations  (3.1-15).  First  a  simplified,  compact  form  of  Equations  (3.1-15)  is  presented. 
This  form  is  useful  when  deriving  the  characteristic  form  of  the  equations.  Because 
Equations  (3.1-15)  are  hyperbolic,  it  is  guaranteed  that  these  equations  are  well-posed  for 
initial  value  problems.  If  these  equations  were  not  well-posed,  any  solution  to  the  initial 
value  problem  would  be  unstable.  This  analysis  is  very  similar  to  the  analysis  performed 
by  Baer  and  Nunziato  [2]  for  their  two-phase  model  equations.  Here  the  same 
characteristic  eigenvalues  are  obtained. 

Though  the  characteristic  form  is  not  immediately  relevant  to  the  work  presented  in  this 
thesis,  it  could  be  important  for  future  work  in  the  unsteady  DDT  problem.  The 
characteristic  form  is  in  some  sense  the  natural  frame  in  which  to  study  the  unsteady 
equations.  The  unsteady  equations  are  transformed  from  a  set  of  partial  differential 
equations  to  a  set  of  ordinary  differential  equations.  Previous  studies  of  the  unsteady 
problem  have  used  the  method  of  lines  to  solve  the  equations  (see  Ref.  47).  With  this 
method  both  time  and  space  derivatives  are  discretized.  Also  to  describe  shock  waves,  it  is 
necessary  to  use  a  special  technique,  such  as  artificial  viscosity  or  flux-corrected-transpon 
(FCT)  to  spread  the  shock  jump  over  a  few  finite  difference  cells.  When  the  characteristic 
form  of  the  equations  is  studied,  no  shock-smearing  method  is  required  to  describe  shock 
jumps. 

This  analysis  will  follow  the  technique  described  by  Whittam  [66]  for  determining  the 
characteristic  eigenvalues  and  eigenvectors.  Consider  a  system  of  partial  differential 
equations  of  the  form 


d\i.  0u. 

A. +  B..'..  ^  =  C. 
Udt  U9x  1 


Multiply  both  sides  of  Equation  (A.l)  by  a  vector  Ij. 


3u.  3u. 

+  l.B..^  =  l.C. 
‘  >J  dt  •  U9x  1  1 


(A.l) 


(A.2) 
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The  vector  Ij  is  chosen  such  that  Equation  (A.2)  can  be  transformed  into  a  system  of 
ordinary  differential  equations.  To  insure  that  Equation  (A.2)  can  be  transformed  to  such  a 
system,  it  is  sufficient  to  require  that  the  following  condition  hold. 

I.B..  =  XIA..  (A.3) 

1  U  1  ij 

where  X,  is  a  variable  scalar  quantity.  If  Equation  (A.3)  holds,  then  Equation  (A.2)  can  be 
written  as 


=  l.C. 
1  1 


(A.4) 


Equation  (A.4)  can  be  transformed  to  an  ordinary  differential  equation  on  special 
curves  in  the  x-t  plane.  On  curves  specified  by 


(A.5) 


Equation  (A.4)  becomes 


du. 

l.A..^  =  l.C. 
»  y  dt 


(A.6) 


To  get  the  form  of  Equation  (A.6)  it  is  required  that  the  eigenvalue  problem  specified  by 
Equation  (A.3)  holds,  that  is 


For  a  non-trivial  solution  to  this  equation  to  exist  it  is  necessary  that 

detf^A,.  -  B..)  =  0 


(A.7) 
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Solution  of  Equation  (A.8)  will  provide  a  set  of  eigenvalues  X.  For  each  eigenvalue,  it 
is  then  possible  to  use  Equation  (A.7)  to  determine  the  vector  Ij.  This  vector  will  have  an 
arbitrary  magnitude.  Using  this  vector  for  the  particular  eigenvalue,  equation  (A.6)  can  be 
used  to  determine  the  characteristic  ordinary  differential  equation  for  the  characteristic 
direction  of  interest  When  substituted  into  Equation  (A.6)  the  arbitrary  magnitude  appears 
as  a  factor  on  both  sides  of  the  equation  and  cancels. 

To  study  the  characteristic  form  of  Equations  (3.1-15),  it  is  first  important  to  write 
these  equations  in  the  reduced  form  required  by  Equation  (A.l).  To  achieve  this  form, 
several  steps  are  necessary.  First,  the  gas  and  solid  mass  equations  are  used  to  eliminate 
density  derivatives  in  gas  and  solid  momentum  and  energy  equations.  Next,  the  reduced 
gas  and  solid  momentum  equations  are  used  to  eliminate  velocity  derivatives  in  the  gas  and 
solid  energy  equations.  Then  the  gas  and  solid  Gibbs  equations  are  used  in  the  gas  and 
solid  energy  equations  to  rewrite  derivatives  of  gas  and  solid  energy  in  terms  of  derivatives 
of  gas  and  solid  entropy  and  density.  Finally,  thermodynamic  relations  developed  in 
Appendix  B  are  used  in  the  gas  and  solid  momentum  equations  to  rewrite  derivatives  of  gas 
and  solid  pressure  in  terms  of  derivatives  of  gas  and  solid  entropy  and  density. 

With  these  steps  and  adopting  Equation  (3.16)  in  favor  of  the  number  conservation 
equation  (3.7),  the  unsteady  two-phase  equations  can  be  written  compactly  as  follows 

(A.9) 


(A.  10) 


a(|ij  3p.  au,  3$.  ap.  4  r  „  i 

3u.  3u.  3(1).  (  8-1  \  9s. 

.<().—  +  p.(().u.-5-^  +  +  <l>.c^-^  +  (J).!  P:  -  = 

i9x  i3x  1  dx  \  *  2  dx 

*  '“2  ■  ’'2'*’l<“2  -  “l>] 


3s.  3s.  3<j).  3<)). 

p.(j).T.^  +  p.(!).T.u.^  -  P.-=^  -  P.u.-^pi.  = 

i9t  1  i9x  1 3t  ‘  ‘3x 


(A.ll) 
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(A.  12) 

3  8r  3  dr  1  1  re, 

-3-  +  -u,^h - — U.-5-4  =  -7C.P. 

r  dt  r  2  3x  dt  Q  2  dx  11 

^2  ^2 

(A.13) 

1  gas  phase  i  =  1  5=1 

with  i 

1  solid  phase  i  =  2  5  =-l 

Equation  (A.9)  represents  the  gas  and  solid  mass  equations;  Equation  (A.  10) 
represents  the  gas  and  solid  momentum  equations;  Equation  (A.  11)  represents  the  gas  and 
solid  energy  equations;  Equation  (A.  12)  is  the  compaction  equation;  and  Equation  (A.  13) 
represents  a  combination  of  the  number  conservation  equation  and  solid  mass  equation. 

The  algebraic  details  required  to  derive  the  characteristic  equations  are  very  lengthy  and 
not  immediately  relevant  to  this  work.  For  this  reason,  only  the  results  will  be  presented 
here.  Six  characteristic  eigenvalues  X  are  found 


X  = 


Uj  ±c 

U2±C 


1 

2 


(A.  14) 


The  characteristics  are  real  and  analogous  to  the  characteristics  found  for  one-phase 
equations. 

The  characteristic  equations  have  been  determined  in  the  limit  when  the  gas  phase  is 
ideal.  There  is  nothing  in  principle  preventing  the  characteristic  equations  from  being 
determined  for  a  non-ideal  gas;  however,  the  algebraic  details  are  much  more  complicated. 
The  characteristic  directions  given  by  Equation  (A.  14)  apply  to  both  ideal  and  non-ideal  gas 
phase  state  equations,  and  the  non-ideal  solid  assumption  has  not  been  relaxed  in  any 
calculations.  Let  Yi  =  and  Yj  =  JC17.  The  equations  in  characteristic  form  in  the  ideal  gas 
limit  are 
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1  ^  ^  1  du.  _ 

1 


5J2-L  p  p’^.^^P/^“2-^)^’'21’,<“2•“l>> 
PiVl  '  ^  ‘  '=i 


+  5' 


<t>. 


Y.p.A.T.r 
‘r  ri  1 


{“lp2P^(Vi■4<“2■^>^■'’/Pi)+“2^(V“P<V^>■"’'3‘^l<V2■■^P'■^) 


5-1  To’' 


It. 


2  8 


2  2 


^2^2 

1  ds.  P.  1  d(b. 

1_L  .  (v-D—i-XjH  = 

V  2  (h  dt.« 

Yi  lO  p.c.  lO 


(A.  15) 


+  5 


0. 


Yp.A.T.r 


1  d<)), 


<f>2  '20 


d<p  r  1  Pi 

=  V,[‘’2■V,-’'l5♦2]-^-i- 

■■  <“20  P.  ‘*'20  ‘  ' 


(A.  16) 
(A.  17) 

(A.18) 


where  the  derivatives  are  defined  as  follows 


i± 


d  _  9  8 

dt.„  8t  ^  8x 
lO 


These  definitions  lead  to  the  following  differential  equations  defining  the  characteristic 
directions 
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—  =  u.  ±c.  on  i±  characteristics 
dt  1  > 


— -  =  u.  on  iO  characteristics 
dt  » 

It  should  be  noted  that  Equations  (A.  15- 16)  reduce  to  familiar  one-phase  formulae 
given  by  Courant  and  Fredrichs  [67]  in  the  one-phase  inert  limit. 
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APPENDIX  B.  Thermodynamic  RELATIONS 


In  this  appendix,  it  will  be  shown  how,  given  a  thermal  equation  of  state  for  pressure 
as  a  function  of  density  and  temperature,  one  can  derive  a  thermodynamically  consistent 
caloric  equation  for  internal  energy  as  a  function  of  density  and  temperature.  This 
technique  will  be  applied  to  the  virial  gas  state  equation  and  solid  Tait  equation.  Equations 
for  sound  speed  and  partial  derivative  of  pressure  with  respect  to  entropy  at  constant 
density  are  derived  for  each  phase.  The  analysis  that  will  follow  is  well-known  in  classical 
thermodynamics  and  can  be  found  in  most  thermodynamics  textbooks. 

General  Analysis 

For  this  analysis  let  the  specific  volume  v  be  defined  as  v  =  1/p.  The  task  is  to  derive 
a  caloric  state  equation  [e  =  e(T,v)]  given  a  thermal  state  equation  [P  =  P(T,v)].  If  energy 
is  to  be  a  function  of  temperature  and  voIudm,  then  the  differential  of  energy  can  be  written 
as  follows: 


de 


dT  + 


de 

8v  T' 


dv 


(B.l) 


The  Gibbs  equation,  Tds  =  de  +  Pdv,  can  be  used  to  write  an  expression  for  the  partial 
derivative  of  energy  with  respect  to  volume: 


(B.2) 


The  specific  heat  at  constant  volume  is  defined  as 


c 

V 


(B.3) 


Equations  (B.2)  and  (B.3)  are  then  substituted  into  Equation  (B.l)  to  yield  the  following: 
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de  =  c  dT  + 

V 


Using  the  Maxwell  relation 


^1  =  ^1 
ar  V  av'r 


(B.4) 


(B.5) 


in  Equation  (B.4)  the  following  equation  is  obtained  for  the  differential  of  energy: 


de  =  c  dT  + 

V 


(B.6) 


which  is  a  convenient  formula  for  determining  a  caloric  equation  given  a  thermal  equation 
of  state. 

Gas  Phase  Analysis 

It  is  assumed  that  the  gas  thermal  equation  of  state  is  given  by 


(B.7) 


By  substituting  Equation  (B.7)  into  Equation  (B.6),  the  following  equation  is  obtained  for 
the  differential  of  gas  internal  energy 


dCj  =  CyjdTj  (B.8) 

By  making  the  assumption  of  a  constant  specific  heat  at  constant  volume,  integrating 
Equation  (B.8),  and  setting  the  arbitrary  integration  constant  to  zero,  the  following  formula 
is  obtained  for  the  gas  internal  energy: 


e  ,  =  c  ,T 

1  vl  1 


(B.9) 
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Internal  energy  can  be  written  in  terms  of  pressure  and  density  by  substituting  Equation 
(B.7)  into  Equation  (B.9)  and  using  the  definition  of  specific  volume. 


"vl 


^  Pj(l+bPj) 


(B.IO) 


The  Gibbs  equation,  TjdSi  =  dej  -  Pi/Pi^  dpj,  can  be  used  with  Equation  (B.IO)  to 
determine  an  expression  for  sound  speed  Cj,  defined  below: 


c 


2 

1 


(B.ll) 


By  using  Equation  (B.IO)  to  determine  the  differential  of  energy  in  terms  of  pressure  and 
density  and  substituting  this  res  ilt  into  the  Gibbs  equation,  the  following  expression  is 
obtained: 


Tjds 


1 


'vl  1 


n - -dP, 

R  p,(i  +  bp,) 


^  p2(l+bpj)2  ' 


(B.12) 


By  holding  entropy  constant  (dsj  =  0),  and  using  Equation  (B.7)  to  reintroduce 
temperature.  Equation  (B.12)  can  be  used  to  determine  an  expression  for  gas  phase  sound 
speed: 


c 


2 

1 


RT  I  1  +  2bp, 


+  (R/c^j)(l 


(B.13) 


It  is  easily  verified  by  setting  b  =  0  that  Equation  (B.12)  reduces  to  the  well-known  ideal 
gas  sound  speed. 

Solid  Phase  Analysis 

For  the  solid  phase  the  assumed  thermal  state  equation  is 
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^2  = 


(Y^-1)Cv^T2 


(B.14) 


By  substituting  Equation  (B.14)  into  Equation  (B.6),  the  following  expression  is  obtained 
for  the  differential  of  solid  energy: 


^^^2  =  ‘'v2'^'^2  (B.15) 

^2 

By  assuming  a  constant  specific  heat  at  constant  volume,  integrating  Equation  (B.15)  and 
assuming  the  arbitrary  integration  constant  is  the  chemical  energy  q,  the  following  equation 
is  obtained; 


®2  "  ‘^v2^2 


P  s 
^20 


V2  +  q 


(B.16) 


Using  the  thermal  state  equation  (B.14)  and  the  definition  of  specific  volume,  Equation 
(B.16)  can  be  rewritten  to  give  internal  energy  as  a  function  of  density  and  pressure. 


e 


2 


‘*2  ^  P,o^ 

(Yj-l)Pj 


+q 


(B.17) 


As  for  the  gas  phase,  the  sound  speed  for  the  Tait  solid  may  be  determined  by 
considering  the  Gibbs  equation.  The  Gibbs  equation  for  the  Tait  solid  in  terms  of 
differential  pressure  and  density,  obtained  from  differentiating  Equation  (B.17),  is 


TjdSj  = 


(Yj  -  1)  Pj 


•dP, 


'’2"V 

(Y,-1)P^ 


(B.18) 


103 


By  setting  the  entropy  change  to  0  in  Equation  (B.18),  solving  for  the  derivative 
representing  sound  speed,  and  using  the  thermal  state  equation  (B.14)  to  reintroduce 
temperature,  the  following  formula  for  the  sound  speed  of  a  Tait  solid  is  obtained: 

4  =  Y,(Yj-1)C^Tj  (B.19) 

Equation  (B.19)  is  identical  to  the  formula  one  finds  for  the  sound  speed  of  an  ideal  gas; 
when  the  sound  speed  is  expressed  as  a  function  of  pressure  and  density,  there  is  a  non¬ 
ideal  term  present. 
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APPENDIX  C.  Model  COMPARISONS:  MOMENTUM  AND  Energy 

EQUATIONS 

In  this  appendix,  the  momentum  and  energy  equations  of  this  work  are  compared  to 
those  of  Baer  and  Nunziato  [2].  The  differences  lie  in  the  particular  form  of  the  pressure 
gradient  term  in  the  momentum  equations  and  in  a  term  known  as  compaction  work  in  the 
energy  equations. 

Momentum  Equations 

The  formulation  of  the  momentum  equations  used  in  this  work,  adopted  from  the  work 
of  Butler  and  Krier  [1],  has  been  criticized  because  it  fails  to  describe  the  equilibrium 
configuration  of  solid  particles  at  rest  in  a  less  dense  fluid  in  the  presence  of  a  gravity  field. 
This  is  not  in  dispute.  It  has  been  argued  that  the  two-phase  equations  as  formulated  by 
Baer  and  Nunziato  are  able  to  describe  such  a  situation  and  thus  are  to  be  preferred  to  the 
Butler-Krier  equations.  Here,  it  is  shown  that  both  formulations  are  in  general  unable  to 
predict  the  equilibrium  situation  described  above. 

The  problem  is  sketched  below  in  Figure  C.l. 


Figure  C.l.  Sketch  of  vertical  settling  problem 


This  sketch  shows  a  mixture  of  a  fluid  and  solid  particles  at  rest  in  a  tube.  Here  a 
gravitational  acceleration,  g,  has  caused  the  heavier  solid  particles  to  settle  to  the  bottom  of 
the  tube. 
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Consider  the  following  two-phase  model  equations,  which  are  inclusive  of  both  the 
Nunziato-Baer  and  Butler-Krier  formulations. 


3uj  9u  9<|) 

=  -8(u,-uj)-p^o^g 


^u. 


au.  ap. 


a<i>.. 


a<i> 


p2'''2'5f*p2V2'^2-3f"'’2-3r^‘'*''’l-5^  =  ■  ®(“2-“l>-p2V 


^ r* 


(C.l) 

(C.2) 

(C.3) 


Here  the  subscript  "1"  represents  the  fluid  phase  and  the  subscript  "2"  represents  the 
solid  phase.  Equations  (C.1)  and  (C.2)  are  the  momentum  equations  for  the  fluid  and  solid 
phases,  respectively.  Equation  (C.3)  is  the  dynamic  compaction  equation.  Density  is 
represented  by  p,  volume  fraction  by  <j>,  velocity  by  u,  pressure  by  P,  drag  coefficient  by  5, 
gravitational  acceleration  by  g,  compaction  viscosity  by  p-g,  and  static  pore  collapse 
function  by  f,  assumed  to  be  a  function  of  only  the  solid  volume  fraction,  <^.  For  k  =  0 
these  equations  describe  the  Nunziato-Baer  formulation,  for  k  =  1,  the  Butler-Krier 
formulation. 

The  following  equations  partially  describe  the  initial  state  in  the  vertical  settling 
problem: 


(j>^(x,0)  =  h(x)  *  (C.4) 

Uj(x,0)  =  U2(x,0)  =  0  (C.5) 


Here,  it  is  assumed  that  there  is  an  initial  distribution  of  particles  given  by  a  general 
function  h(x).  It  is  further  assumed  that  both  particles  and  fluid  are  at  rest. 

It  would  seem  that  a  basic  test  for  any  model  of  this  problem  is  that  the  model  should 
predict  that  the  mixture  stays  at  rest;  thus  at  this  initial  state,  the  equations  should  predict 
that  no  variables  change  with  respect  to  time.  To  insure  that  no  volume  changes  are 
predicted,  a  condition  on  the  relation  between  P2,  Pi  and  ((>2  is  obtained  from  Equation 

(C.3): 


=  Pj+f[h(x)] 
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(C.6) 


An  equilibrium  condition  is  also  obtained  from  the  fluid  momentum  equation  (C.  1)  by 
using  the  initial  conditions  (C.4,  5),  For  no  fluid  motion  to  be  predicted,  the  following 
condition  must  hold: 


3P,  Pi  dh 

ST  “  ■  PiS  <C.7) 

For  K  =  0,  a  result  identical  from  one-phase  fluid  statics  is  recovered.  It  has  been 
argued  [68]  that  this  limit  must  also  be  recovered  from  a  two-phase  model  and  that  this  is  a 
sufficient  reason  to  take  k  =  0.  However,  it  is  still  not  clear  whether  this  familiar  result 
should  extend  to  the  two-phase  situation. 

When  the  solid  momentum  equation  (C.2)  is  considered,  it  is  seen  that  both 
formulations  have  difficulty  describing  an  equilibrium  configuration.  By  substituting  the 
initial  conditions  (C.4,  5),  the  compaction  equation  condition  (C.6),  and  the  fluid  pressure 
gradient  condition  (C.7)  into  Equation  (C.2),  the  following  condition  is  obtained  for  no 
solid  acceleration: 


(p/P2-l)g-J-^hf(hW)].  ic 


1  1 


P  h 


p  h  dx 


^.0 


(C.8) 


This  equation  raises  questions  regarding  equilibrium  in  the  presence  of  gravitational 
forces  and  initial  volume  fraction  gradients..  Consider  two  limiting  cases  for  the  Nunziato- 
Baer  model  (k  =  0),  In  the  first  case  consider  a  situation  in  which  there  is  no  intragranular 
stress;  the  particles  are  in  contact  but  are  not  exerting  a  force  on  each  other.  This  would 
correspond  to  the  condition  f  =  0.  In  this  limit,  equation  (C.8)  predicts  equilibrium  only 
when  the  fluid  density  is  equal  to  the  solid  density,  which  in  general  is  not  satisfied.  (In 
this  case  it  is  questionable  whether  the  state  equations  would  allow  Equation  (C.6)  to  hold 
also.)  In  the  second  case  consider  the  zero-gravity  limit,  g  -» 0.  In  this  limit  for  Equation 
(C.8)  to  be  satisfied,  the  static  pore  collapse  function  f  must  be  of  the  form  f  =  constant  / 
(1)2,  a  condition  which  is  not  enforced  by  the  Nunziato-Baer  model.  The  condition  (C.8) 
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has  not  been  enforced  by  the  Butler-Krier  model.  It  should  be  said  that  the  Butler-Krier 
model  never  attempts  to  describe  the  situation  in  Figure  C.  1. 

From  this  analysis,  it  is  clear  that  both  the  Nunziato-Baer  and  Butler-Krier  model 
equations  are  incapable  of  describing  an  equilibrium  state  in  the  presence  of  either  a  volume 
fraction  gradient  or  gravity  forces  unless  very  restrictive  conditions  are  placed  on  the 
constitutive  relations.  Neither  of  these  models  currendy  enforces  such  restrictions. 

Energy  Equations 

Baer  and  Nunziato  have  included  a  term  in  their  model  which  is  intended  to  model 
experimentally-observed  hot  spots  in  granular  explosives  and  the  work  associated  with  the 
local  distortion  of  grains  when  a  granular  material  is  compacted.  This  term,  called 
compaction  work,  appears  in  both  the  solid  and  gas  phase  energy  equations.  It  is 
constructed  such  that  when  compaction  work  is  predicted,  energy  is  removed  from  the 
solid  phase  and  deposited  in  the  gas  phase.  This  local  energy  deposition  gives  rise  to  a 
local  hot  spot  which  encourages  a  local  acceleration  of  the  reaction  rate.  It  is  shown  by 
Baer  and  Nunziato  that  this  compaction  work  term  is  consistent  with  but  not  required  by  the 
second  law  of  thermodynamics. 

Here,  it  will  be  shown  that  the  presence  of  compaction  work  gives  rise  to  a 
fundamental  inconsistency  in  the  limit  of  an  inert  mixture  where  the  ratio  of  initial  gas 
density  to  initial  solid  density  is  small.  In  this  limit  the  steady-state  mixture  energy 
equation  predicts  a  result  inconsistent  with  the  solid  energy  equation.  It  is  shown  that  the 
solid  energy  equation  in  this  limit  gives  rise  to  energy  escaping  from  the  system. 

Consider  the  following  equations,  general  equations  which  encompass  the  gas  and 
solid  energy  equations  of  both  models. 

P,*.(e,+u?/2)J  .  8P,u,^  = 

■  h(Tj-Tj)  +  o(Uj-Uj)Uj  -  +  5^P2-‘’r«VX'’2'^“'’2>) 

^ r 


(C.9) 
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Pj<l>i‘2«^)]  +4[p2'l>2“2{'2*'4oK«j“2l  +  SP,»; 


a<l>i 

2  ax 


h(Tj-T2)  -  aCuj-u^Uj  +  V“^)  - 


10) 


5  = 


'  0  Model  of  this  Work 

1  BaCT-Nunziato  Model 


The  new  notation  introduced  here  is  that  h  is  considered  to  be  a  general  function 
specifying  the  heat  transfer  coefficient,  likewise  a  is  taken  to  be  a  general  function 
specifying  the  drag  coefficient,  and  c^^  is  a  general  function  specifying  the  combustion  rate. 
The  parameter  5  is  used  to  distinguish  the  two  model  formulations. 

When  Equations  (C.IO)  and  (C.ll)  are  added,  a  homogeneous,  unsteady  mixture 
energy  equation  is  obtained. 

|(p,'i>i(=i+“?«)+Pjfj(y"H]  * 

It  is  argued  by  Baer  [45]  that  in  the  limit  where  material  is  inert  =  0)  and  where  the 
effect  of  the  gas  phase  is  negligible  that  Equation  (C.IO)  reduces  to  the  following 

p2l'i'2-^“^)]+|[p2'!'2“2('2«^)*P2*2"2]  = 

(C.12) 

If  Equations  (C.  1 1)  and  (C.  12)  are  transformed  to  steady  dimensionless  form  and  the  limit 
of  zero  gas  phase  density  is  taken,  the  inconsistency  becomes  apparent.  Using  the  same 
technique  and  nomenclature  found  in  the  main  text  for  writing  steady  dimensionless 
equations,  it  is  found  that  Equation  (C.ll)  transforms  to  the  following  (equivalent  to 
Equation  (5.7)) 
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ll5P,<l>,v,l 


P 


■^PjVa 


2  p 
'^2  ^2 

^2/ 


Vll^^6^14  +  1/2  +  TCj^)  +  (1  -  7Cjj)(7tj4  +  TCg  +  TCjq  +  1/2  +  7C2j)  (C.13) 


In  the  limit  of  zero  gas  phase  density,  JC5  is  zero.  Thus  in  this  limit  the  steady  mixture 
equation  (C.13)  becomes 


PaVa 


2  p 

^2"TV 

^2/ 


(l-7lii)(7Ci4  +  7C8  +  7tio+l/2  +  7l2i) 


(C.14) 


The  steady  dimensionless  form  of  Equation  (C.12)  is 


_d_ 


P2V2 

. 

,^2*Tr'’"V2 

=  5lt,(|)|(|>j(P2-«i|'j))  (C.15) 


It  is  obvious  that  Equations  (C.14)  and  (C.15)  are  consistent  only  when  5  =  0,  that  is  when 
compaction  work  is  ignored. 

Inclusion  of  compaction  work  leads  to  violation  of  the  conservation  of  energy  in  the 
zero-density  gas  phase  limit.  This  is  easily  seen  by  considering  the  application  of  the 
unsteady  energy  equation  (C.12)  to  the  following  problem  (see  Fig.  C.2).  Strike  with  a 
piston  a  constant  area  tube  closed  at  one  end  containing  a  porous  material.  After  a  period  of 
time  bring  the  piston  to  rest.  The  piston  motion  induces  a  pressure  imbalance  in  the  porous 
material  (i.  e.  (P  -  f)  >  0).  After  the  piston  is  brought  to  rest,  a  zero  velocity  boundary 
condition  must  be  enforced  at  both  ends  of  the  tube.  However  the  material  inside  the  tube 
is  not  in  a  state  of  equilibrium. 
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Figure  C2  Sketch  of  Piston  Problem 

By  integrating  the  energy  equation  (C.12)  from  x  =  0  to  x  =  L,  it  is  seen  that  for  5  =  0,  the 
energy  of  the  system  is  conserved  and  for  8  =  1,  energy  leaves  the  system  as  time 
progresses.  The  time  rate  of  change  of  energy  per  unit  cross-sectional  area  for  this  system 
is 


0 

For  0  <  <1)2  ^  1,  the  integrand  of  the  right  hand  side  of  the  total  energy  equation  is  always 
positive;  therefore,  for  8  =  1,  energy  leaves  this  system,  and  for  5  =  0,  energy  is 
conserved.  In  order  to  preserve  energy  conservation  in  this  limit,  and  in  light  of  the  fact 
that  compaction  work  is  not  necessarily  required  by  the  second  law  of  thermodynamics, 
compaction  work  is  not  included  in  this  model. 

It  is  concluded  that  though  it  may  be  important  to  model  hot  spot  formation,  the 
proposed  mechanism  of  compaction  work  has  an  inherent  flaw,  and  in  order  to  model  such 
phenomena  another  model  must  be  proposed.  To  model  hot  spots  in  a  granular  material 
which  arise  from  the  material  compaction  is  difficult  in  the  context  of  a  two-phase  mixture 
model.  One  would  need  to  devise  a  way  to  non-uniformly  distribute  the  energy  introduced 
to  the  system  by  the  piston  (P  dV  work)  to  the  particles.  The  non-uniform  distribution 
would  allow  some  particles  to  have  higher  temperatures  than  others,  thus  giving  rise  to 
local  "hot  spots."  It  is  unclear  how  this  can  be  achieved  with  a  two-phase  mixture  model 
which  relies  on  averaged  properties.  In  fact  one  of  the  strengths  of  two-phase  modeling  is 
that  details  of  microstructure  do  not  need  to  be  considered  as  these  local  variations  are 
eliminated  in  the  averaging  process.  For  this  reason,  it  may  be  impossible  to  attempt  to 
describe  hot  spots  with  a  two-phase  mixture  model. 
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Appendix  d.  two-Phase  CJ  Deflagrations 

It  is  possible  in  principle  to  use  the  two-phase  model  to  study  two-phase  deflagrations, 
reactive  waves  which  travel  much  slower  than  detonations  and  which  have  a  much  lower 
pressure,  temperature,  and  density  rise.  Because  of  the  more  moderate  changes  in  the  state 
of  the  gas,  the  ideal  gas  state  equation  is  appropriate  for  use  in  studying  two-phase 
deflagrations.  Understanding  of  two-phase  deflagrations  can  be  gained  by  studying  the 
complete  reaction  two-phase  Rayleigh  line  and  Hugoniot  equations  (5.21,  22)  in  the  Pj- 
1/pj  plane  (see  Fig.  D.l). 


Figure  D.l  Two-Phase  Complete  Reaction  Deflagration  Rayleigh  Line  and  Hugoniot 

Deflagration  solutions  are  found  at  the  intersection  of  these  two  curves  at  gas  pressures 
lower  than  the  initial  apparent  pressure  P^  and  gas  densities  lower  than  the  initial  apparent 
density  p^.  It  is  possible  to  predict  a  maximum  deflagration  wave  speed,  called  here  the  CJ 
deflagration  speed.  At  the  CJ  deflagration  state,  the  two-phase  Rayleigh  line  is  tangent  to 
the  two-phase  Hugoniot.  For  wave  speeds  greater  than  the  CJ  deflagration  speed,  but  less 
than  the  CJ  detonation  speed,  there  is  no  intersection  of  the  two-phase  Rayleigh  line  and 
Hugoniot  and  thus  no  solution.  For  wave  speeds  less  than  the  CJ  deflagration  speed  two 
solutions  are  obtained,  a  strong  and  weak  deflagration  solution. 

For  an  ideal  gas  the  complete  reaction  CJ  wave  speed  is  given  exactly  by  the  following 
equation  using  the  nomenclature  of  Chapter  5. 
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(2c^j+R)  -  cJjPyp^  ±  J (2c,j+R{ey(2c^,+R)  -  (2e^+P/p^)cJ,P^R/pJ 

°CJ  i 

(D.l) 

Here  the  plus  branch  of  this  equation  corresponds  to  the  CJ  detonation  state  and  the 
minus  branch  corresponds  to  the  CJ  deflagration  state.  When  Pa/(Pa®a)  «  ^he  CJ 
deflagration  state  simplifies  considerably.  In  this  limit,  which  is  relevant  for  many  physical 
systems  including  the  HMX  system  studied  in  this  thesis,  the  CJ  deflagration  state  is 
approximated  by  the  following  equations. 
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(D.3) 

(D.4) 

(D.5) 

(D.6) 

(D.7) 


It  is  important  to  stress  that  beyond  describing  the  maximum  speed  two-phase 
deflagration  wave  the  interpretation  of  Equations  (D.2-5)  is  unclear.  At  this  point  it  is  not 
known  whether  a  steady  two-phase  deflagration  structure  can  be  predicted  by  the  model 
equations  (5.1-15)  and  if  such  a  structure  could  be  predicted,  what  conditions  would  dictate 
whether  a  CJ,  strong,  or  weak  deflagration  was  obtained.  A  limited  study  was  undertaken 
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to  find  steady  deflagration  structure  with  no  success.  This  study  included  CJ  deflagrations 
along  with  strong  and  weak  deflagrations.  Kuo  and  Summerfield  have  found  steady  two- 
phase  deflagration  structvire  using  a  similar  model  [37].  Also  both  the  Kuo  and 
Summerfield  model  and  the  model  of  this  work  have  neglected  diffusive  processes  such  as 
heat  conduction  and  viscous  momentum  transport  which  may  be  very  important  for  the 
relatively  slow  deflagration  waves. 
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APPENDIX  E.  DERIVATION  OF  UNCOUPLED  EQUATIONS 


This  appendix  will  provide  a  detailed  explanation  of  how  the  coupled  set  of 
differential-algebraic  equations  describing  steady  detonation  structure  (5.1-15)  can  be 
written  as  four  differential  equations  in  four  unknowns  and  how  in  the  zero  heat  transfer, 
zero  compaction  limits  these  equations  can  be  further  reduced  to  form  the  two-equation 
model  (5.45-46).  First,  it  will  be  shown  how  through  an  algebraic  analysis  the  mixture 
equations  (5.5-7)  can  be  used  to  write  gas  phase  quantities  in  terms  of  solid  phase 
quantities.  It  is  found  that  this  process  involves  the  solution  of  a  cubic  equation.  Next  the 
coupled  differential  equations  (5.1-4)  are  uncoupled  using  linear  algebra  techniques.  It  is 
then  shown  how  these  equations  reduce  to  the  two-equation  model. 

The  mixture  equations  (5.5-7),  solid  and  gas  caloric  state  equations  (5.13,  10),  and 
porosity  definition  (5.15)  are  rewritten  here 
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By  using  Equations  (E.4)  and  (E.5)  to  eliminate  gas  and  solid  energy  in  Equation 
(E.3)  and  Equation  (E.6)  to  eliminate  gas  volume  fraction  in  Equations  (E.1-3),  Equations 
(E.l -3)  can  be  rewritten  as  follows 
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P.''i=a(p/V''2-'’2) 

P,+PvJ  =  b(p2.'I>2.V2.P2) 
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where  A,  B,  C,  and  D  are  functions  of  solid  density,  volume  fraction,  velocity,  and 
pressure  defined  below 
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Equations  (E.7-9)  can  be  combined  to  form  a  cubic  equation  for  gas  density.  This  is 
done  by  first  using  Equation  (E.7)  to  express  gas  velocity  as  a  function  of  gas  density  and 
solid  variables.  Then  gas  velocity  may  be  eliminated  from  Equations  (E.8)  and  (E.9). 
Thus  modified.  Equation  (E.8)  can  be  used  to  express  gas  pressure  as  a  function  of  gas 
density  and  solid  variables.  This  result  is  used  to  eliminate  gas  pressure  from  the  modified 
energy  equation  (E.9).  The  resulting  equation  is  a  cubic  equation  for  gas  density  whose 
coefficients  are  functions  of  the  solid  phase  density,  volume  fraction,  velocity,  and 
pressure. 
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2(7t^  -  l)7Cj3Cj  p3  +  -  1)(ABiCj3  -  C) j  p2  + 

[^2AB7i.^-(7i^-l)A^7cJp^  -[(Tt^+DA^j  =  0  (E.12) 

Equation  (E.12)  can  be  solved  exactly  for  gas  density  in  terms  of  solid  phase  variables 
and  parameters.  The  solution  is  very  lengthy  and  can  be  easily  produced  using  the  formula 
for  solution  to  the  general  cubic  equation.  Three  roots  are  found  for  Equation  (E.12).  One 
is  associated  with  a  shocked  gas  and  is  analogous  to  the  strong  branch  of  the  ZND 
solution.  Another  is  associated  with  an  unshocked  gas  and  is  analogous  to  the  weak 
branch  of  the  ZND  solution.  The  third  predicts  a  negative  density  for  all  cases  studied  and 
is  rejected  as  unphysical.  This  root  is  not  present  when  non-ideal  gas  effects  are  absent. 
(It  is  seen  from  Equation  (E.12)  that  for  no  non-ideal  effects,  7ti3  =  0,  that  the  equation  is 
quadratic,  and  only  two  roots  are  present.)  It  is  possible  for  Equation  (E.12)  to  predict  a 
pair  of  imaginary  roots  under  certain  conditions.  If  such  a  condition  was  reached,  the 
detonation  structure  must  be  rejected  as  unphysical.  In  addition  to  solving  for  the  gas 
phase  variables  within  the  reaction  zone.  Equation  (E.12)  is  used  to  determine  the  shock 
state  of  the  gas. 

With  the  gas  density  predicted  from  Equation  (E.12)  as  a  function  of  solid  phase 
variables,  all  other  gas  phase  variables  can  be  expressed  as  functions  of  solid  phase 
variables.  The  gas  velocity  is  found  by  using  Equation  (E.7).  The  gas  pressure  can  then 
be  determined  from  Equation  (E.8)  and  the  energy  from  the  state  equation  (E.5).  The  gas 
temperature  and  sound  speed  can  then  be  found  using  Equations  (5.9,1 1). 

In  the  numerical  code  which  predicts  reaction  zone  structure.  Equation  (E.12)  was 
solved  using  the  IMSL  subroutine  ZRPOLY.  Though  one  could  use  the  exact  cubic 
solution  to  determine  the  gas  density,  the  numerical  accuracy  of  the  solution  is  higher  when 
ZRPOLY  is  used.  Given  a  general  polynomial  equation,  the  subroutine  ZRPOLY 
determines  all  roots,  real  and  complex. 

Equations  5. 1-4  can  be  expressed  in  the  form 

du. 

A..^  =  B.  (E.13) 

‘Jd^ 

where  Uj  =  (p2,  <|)2,  V2,  P2)  and  Ajj  and  Bj  are  functions  of  P2,  (1)2,  ^2,  and  P2.  To  put  the 
equations  in  a  form  suitable  for  phase  space  analysis,  explicit  expressions  for  the 


I 
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deri\  atives  dUj/d^  must  be  obtained.  This  is  done  by  multiplying  both  sides  of  Equation 
(E.13)  by  the  inverse  of  Ay. 


du.  , 

-1  =  a:.'  b. 


One  necessary  step  to  express  Equations  (5.1-4)  in  this  form  is  to  use  the  solid  caloric 
state  equation  (5.13)  to  determine  an  expression  of  the  derivative  of  solid  energy  in  terms 
of  solid  pressure  and  density.  This  derivative  is  given  below: 


^  ^  1  ^2  .  P2  ^17^8  % 

d^  (jCj^  - 1)  d^  (Jt^^  -  1)  p2  d^ 


(E.14) 


Expanding  the  derivatives  in  Equations  (5.1-4)  and  using  Equation  (E,14)  allows  a 
system  in  the  form  of  Equation  (E.13)  to  be  written. 
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The  left  side  of  Equation  (E.15)  is  expressed  in  terms  of  the  fundamental  variables  P2,  <))2, 
V2,  and  P2.  The  right  hand  side  can  also  be  expressed  in  terms  of  these  variables.  The 
method  described  earlier  in  this  appendix  can  be  used  to  write  the  gas  phase  variables  Vj, 
Pj,  and  Tj  as  functions  of  the  fundamental  variables.  The  number  conservation  equation 


(5.8)  and  solid  thermal  state  equation  (5.12)  can  be  used  to  express  the  radius  r  and  solid 
temperature  T2  as  a  functions  of  the  fundamental  variables. 

By  multiplying  both  sides  of  Equation  (E.15)  by  explicit  expressions  can  be 
obtained  for  the  derivatives  of  the  fundamental  variables. 
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where  D,  E,  and  F  are  defined  as  follows 
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F  =  -;c3(ic,T2- 


(E.22)  j 
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Equations  (E.16, 18, 19)  are  singular  when  the  following  condition  is  met 
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By  using  the  solid  thermal  state  equation  (5.12)  to  eliminate  solid  pressure  and  density 
in  favor  of  solid  temperature  and  using  the  solid  sound  speed  definition  (5.14)  it  is  seen 
that  Equation  (E.23)  can  be  rewritten  as 
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Thus  when  the  velocity  of  the  solid  relative  to  the  wave  head  is  locally  sonic,  the  system  of 
equations  (E,15)  is  singular.  It  is  seen  by  examination  of  Equations  (E.21-22)  that  the 
equations  are  also  singular  at  the  complete  reaction  state  because  the  particle  radius  r 
appears  in  the  denominator  of  the  expressions  for  E  and  F. 

The  two-equation  model  can  be  derived  form  Equations  (E.  16-22).  To  derive  the  two- 
equation  model,  one  must  consider  the  zero-compaction,  zero-heat  transfer  limit, 
corresponding  to  %  0, 7C3  0.  In  this  limit.  Equations  (E.20,  22)  hold  that  D  =  F  =  0. 

Then  if  Equation  (E.16)  is  multiplied  by  solid  velocity  and  added  to  the  product  of  solid 
density  and  Equation  (E.  18),  the  following  homogeneous  equation  is  obtained. 

dp  dv. 

This  equation  can  be  integrated  to  form  the  algebraic  relation 


PjVj  =  ■ ' 


(E.26) 
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In  applying  the  initial  conditions  when  integrating  this  equation,  it  is  unnecessary  to 
specify  whether  the  initial  state  is  shocked  or  unshocked.  This  arises  because  Equation 
(E.26)  is  equivalent  to  the  shock  relation  (5.29)  when  it  is  considered  that  porosity  does  not 
jump  through  the  shock  wave. 

To  determine  a  second  algebraic  relation.  Equation  (E.26)  must  first  be  used  to 
eliminate  solid  velocity  in  favor  of  solid  density  in  all  remaining  equations.  Then  if 
Equation  (E.19)  is  multiplied  by  the  factor 


1 


and  added  to  the  product  of  Equation  (E.16)  and  the  factor 
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the  following  homogeneous  equation  is  obtained 
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This  equation  may  be  integrated  to  form  an  algebraic  relation  between  solid  pressure 
and  density.  The  constant  of  integration  for  this  expression  is  dependent  on  whether  the 
initial  state  is  shocked  or  unshocked. 
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When  Equations  (E.26)  and  (E.28)  are  used  to  eliminate  solid  velocity  and  pressure 
from  Equations  (E.16)  and  (E.17),  the  two-equation  model  is  found. 
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By  multiplying  the  numerator  and  denominator  of  the  right  side  of  Equation  (E.30)  by 
the  factor 
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the  form  of  Equations  (5.45, 46)  is  found. 
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APPENDIX  F.  DERIVATION  OF  NUMBER  CONSERVATION  EQUATION 


As  the  number  conservation  is  not  universally  used  in  two-phase  granular  detonation 
theory,  a  derivation  for  the  number  conservation  equation  (3.7)  along  with  Equation  (3.16) 
is  given  here. 

The  volume  fraction  of  particles  <1)2  is  defined  as  the  ratio  of  the  volume  of  particles  to 
the  total  volume. 


Volume  Particles 
^2  ~  Total  Volume 


(F.l) 


If  it  is  assumed  that  the  particles  are  spheres  of  radius  r,  then  the  volume  of  particles  is 
equal  to  the  product  of  the  number  of  particles  and  the  volume  of  a  single  particle.  Based 
on  this  assumption  Equation  (F.l)  is  written  as 


<t>  = 


4  3 

(Number  of  Particles)  —  nr 
_ 3 _ 

Total  Volume 


(F.2) 


If  the  number  density  n  is  defined  as  the  number  of  particles  per  total  volume,  then 
Equation  (F.2)  can  be  used  to  determine  an  expression  for  number  density  as  a  function  of 
particle  radius  and  solid  volume  fraction. 


(F.3) 


By  performing  a  simple  control  volume  analysis,  an  expression  for  the  conservation  of 
number  density  can  be  derived.  The  expression  is 
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This  equation  could  be  modified  to  describe  the  agglomeration  or  splitting  of  particles  by 
use  replacing  the  zero  on  the  right  side  of  Equation  (F.4)  with  a  functional  relation  suitable 
to  describe  such  a  phenomenon. 

Substituting  the  number  density  definition  (F.3)  into  the  number  conservation  equation 
(F.4),  an  equation  identical  to  Equation  (3.7)  is  derived. 

=  0  <F-5) 


Using  the  Galilean  transformation  ^  =  x  -  Dt,  V2  =  U2  -  D  where  D  is  the  steady  wave 
speed  allows  Equation  (F.5)  to  be  written  in  steady  form. 


=  0 


(F.6) 


Using  the  initial  conditions  from  Chapter  3,  Equation  (F.6)  may  be  integrated  to 
provide  the  following  algebraic  expression  for  particle  radius  as  a  function  of  particle 
velocity  and  volume  fraction 


(F.7) 


To  obtain  an  explicit  equation  for  the  particle  radius  evolution.  Equation  (F.5)  can  be 
expanded. 


(F.8) 


The  solid  mass  equation  (3.2)  can  be  used  to  write  an  expression  for  the  derivative  of 
solid  volume  fraction. 
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By  using  Equation  (F.9)  to  eliminate  the  volume  fraction  derivative  in  Equation  (F.8), 
an  expression  identical  to  Equation  (3.16)  fcff  the  evolution  of  particle  radius  is  obtained. 


ar  ar 


(F.IO) 


This  equation  states  that  the  particle  radius  changes  in  response  to  combustion  and 
compressibility  effects.  Equation  (F.IO)  is  inconsistent  with  the  model  equation  used 
formerly  by  Krier  and  co-woricers  to  determine  the  particle  radius.  As  stated  in  Ref.  1,  the 
particle  bum  law  used  in  these  works  is  (correcting  for  a  sign  error  in  the  paper) 


dr 
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(F.ll) 


In  this  equation  the  definition  of  the  derivative  d/dt  is  unclear  as  to  whether  or  not 
convective  terms  are  included.  Regardless  of  this  question,  it  is  clear  that  Equation  (F.ll) 
does  not  account  for  compressibility  effects  in  the  particles.  It  must  be  concluded  that  since 
Equation  (F.ll)  is  inconsistent  with  Equation  (F.IO)  that  the  model  of  Ref.  1  does  not 
conserve  number,  thus,  the  physical  motivation  of  Equation  (F.l  1)  is  unclear. 


